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ABSTRACT 

The present work is a Kalman filter study, in indirect 
feedback configuration, for a proposed integrated inexpen- 
sive Inertial Navigation System/Global Positioning System 
(I.N.S./G. P.S. ) . 

A one nautical mile per hour, local-level, two- 
accelerometer I.N.S. is used where the errors are 
represented by a 7 state linear model . 

G.P.S. is assumed to provide four range measurements 
from an equal number of satellites with the best relative 
position among those in view. 

I.N.S. error analysis showed error dependence on 
Schuler frequency and that it was possible to neglect 
Foucault modulation for navigation purposes. 

The present I . N . S . /G . P . S. system has been shown to be 
quite effective since the navigation errors are reduced 
quickly for both short and long term periods without any 
divergence . 
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I. INTRODUCTION 



A. OVERVIEW OF AN INERTIAL NAVIGATION SYSTEM (I.N.S.) 

A conventional gimballed inertial measurement unit 
consists of a platform suspended by a gimbal structure that 
allows three degrees of rotational freedom [Ref. 1,2,7]- 
The outermost gimbal can be attached to the body of some 
vehicle and allow that vehicle to undergo any change in 
angular orientation while maintaining the platform fixed 
with respect to some desired coordinate frame. 

Gyros mounted on the platform sense the angular rate of 
the platform with respect to inertial space and their 
outputs are sent through electronics to the torque motors on 
the gimballed structure, commanding them to maintain a 
desired platform orientation regardless of the orientation 
of the outermost gimbal which remains fixed to the body. 

Feedback control loops that keep the gyro outputs 
nulled , will maintain at the same time the platform fixed 
with respect to the inertial space. These feedback loops 
are such that, in practice, the platform orientation is kept 
essentially stable regardless of the most violent vehicle 
maneuvering. Additional (computed) inputs can be added to 
the above feedback loops to maintain some other orientation. 
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such as North-East-Down, correspond ing to the current 
location of the vehicle. 

Accelerometers mounted on the platform can provide the 
vehicle's acceleration with respect to the known set of 
reference coordinates. In fact, specific force is measured 
by the accelerometers so that local gravity must be computed 
and appropriately subtracted from these sensor outputs 
in order to obtain a measurement of actual vehicle 
acceleration . 

The vehicle's velocity and position are obtained by 
integration of the above acceleration measurement signals. 
Attitude information as well as translational information is 
provided by the I.N.S.. A typical gimballed inertial 
measurement unit [Ref. 2] is shown in Fig. 1. 

B. OVERVIEW OF THE G.P.S. 

The Global Positioning System (GPS) is a satellite 
navigation system currently under development. It will 
consist, according to today's available information, of 18 
satellites placed in groups of six in each of three 

different circular, 12 hour orbits at an altitude of 10,900 

o o 

N.M. inclined 63 to the equator and spaced 120 apart. 

The satellites will broadcast pseudo-random noise codes 
(codes P and C/A) and ephemerides on two L-band signals to 
users worldwide in such a way that each satellite signal can 
be distinguished from the others by the user. A user will 
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Gimbal two 




Figure 1. Typical Gimbaled Inertial Measurement Unit [Ref. 2] 
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be equipped with a small receiver (G.P.3. user equipment) 
which measures the pseudo-range and pseudo-range rate from 
the user to the satellite . 

By means of a correlator -detec tor the time (phase) shift 
between each satellite signal and the user's unsynchronized 
clock will be measured in his receiver to provide an 
indication of the range from the satellite to the user. 
Typically, four satellite signals may be received simul- 
taneously by the user equipment. 

The phases of the NA VSTAR/G . P. S. system are shown in 
Fig . 2 C Ref . 8 ] . 

C. I.N.S. OPTIMAL AIDING 

Once we have available a typical inertial measurement 



unit, or 


the i 


ner t 


ial nav 


igation sys 


tern as 


a 


whole , 


t 


he 


question 


natur 


all y 


ar ises 


: why does 


this 


system re 
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ires 


optimal 


aiding 


by 


other n 
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> The 


an 


swer to 


this que 


stion 


is g 


iven in 


Ref. 2 and 


here 


we 


pr e se n 


t 


the 


concepts 


only . 



















Due to the tight control loops supporting the I.N.S. 
very good high frequency information is provided. However, 
because of gyro characteristics, the system drifts at a slow 
rate so that the long term (low frequency content) of the 
data is poor . It is well known that all inertial systems 
have position errors that grow slowly with time and these 
errors are unbounded. 



22 




00 



c/5 

O 

c/i 

CS 

E 



CO 



o 



< 

CO 

> 

< 



CNI 

o 

^1 

oo 

• H 

P-. 



23 



As opposed to an I.N.S. which can be classed as a "one 



nautical mile per hour system” due to the associated 
position error, most other navigation aids provide very good 
low frequency information but subject to considerable high 
frequency noise, due to instrument noise, atmospheric 
effects, antenna oscillation, unlevel ground effects and so 
forth . 

One would want to combine the available information from 
an I.N.S. and other external sources in an optimal manner if 
possible so that one can obtain efficient estimates of 
navigation parameters that are best with respect to some 
well defined criterion. Such an optimal approach is 
provided by the Kalman filter approach which is briefly 
discussed next. 

D. KALMAN FILTER 

The Kalman filter is an optimal recursive data process- 
ing algorithm located in the on-board computer or central 
processor that uses sampled data with sample period on the 
order of 5-60 seconds, to maintain estimates of approxi- 
mately 60-70 state variables. The filter combines all 
available measurement data with prior knowledge of the 
system and measuring devices to produce an estimate of the 
system states in such a manner as to statistically minimiie 
the resulting errors. In more easily understood terms the 
filter, or computer program, uses the statistical 
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characteristics of the errors in both the inertial naviga- 
tion components and the external information providing the 
best estimate possible, subject to certain modeling 
assumptions • 

The filter will act to optimize the attitude, position, 
and velocity information accuracy by weighting each data 
source heavily in the frequency ranges where it provides 
more accurate information, and suppressing it in the region 
where it is less accurate. The inertial system provides 
good high frequency information but it drifts slowly and 
therefore exhibits poor low frequency performance. On the 
other hand, the external aids (such as G.P.S.) generally 
exhibit good low frequency information but are subject to 
high frequency noise. Therefore, the filter will use the 
good low frequency external (G.P.S.) information to damp ou 
the slowly growing errors in the inertial system. 

1 . Type of Filter Implementation 

There are two very important aspects of implemen- 
tation of a Kalman filter in conjunction with an inertial 
system [ Ref . 2 ] . 

a) Total state space (direct) versus error state 
space (indirect) formulation, and 

b) Feedforward versus feedback mechanizations. 

In the indirect formulation the errors in the I.N.S. 
indicated position and velocity are among the estimated 
variables and each measurement presented to the filter is 
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the difference between the I.N.S. and the external source 
(G.P.S.) data. The I.N.S. itself follows the high frequency 
motions of the vehicle very accurately, and there is no need 
to model these dynamics explicitly in the filter but the 



dynamics upon win 


lich the filter is 


based is a 


set 


of 


inertia 


system 


err 


or propagation 


equations 


, which are 


r elat 


ivel y 


well d 


evel 


oped , 
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ved , low 


frequency , 
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very 
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linear [ 


Ref. 2, pp. 


296 ] . 
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adequacy of 
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lin 


ear mo<je 



is enhanced. Second is the fact that many of the predicted 
error states which at next time sample time are zero, need 
not be computed explicitly. 

The indirect feedback configuration of the Kalman 
filter is shown in Figure 3* In Ref. 2 there is explicitly 
documented the discussion of the Kalman filter configura- 
tions and mechanizations with their advantages and more 
comments . 

2. Assumptions 

The Kalman filter can be shown to be the best filter 
of any possible form based on the following three 
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Figure 3. Indirect Feedback Kalman Filter |Ref. 2] 




assumptions: a linear system, white noise drivers, and the 



Gaussian distribution of noise . 

Although the system itself may be nonlinear, 
formulation of an approximate linear error state space model 
makes linear analysis possible. The justification for the 
linear model is based on two points. For the aided I.N.S. 
case the use of linear error state space models yields a 
very adequate representation. The techniques of linear 
system analysis are also well developed and better 
understood than those of nonlinear analysis. 

The white noise assumption implies that the noise 
is not correlated in time and also has equal power at all 
frequencies. If, in fact, a time correlated noise is 
required to adequately model the system, it can be produced 
oy passing white noise through a linear shaping filter. The 
system can then be modeled with an augmented state variable 
as a linear system driven by white noise. 

Gaussianess pertains to the distribution of 
amplitudes of the noise and implies that at any single point 
in time the probability density function of the amplitude 
takes on the shape of a normal bell-shaped curve. The 
assumption of Gaussian noise amplitude is justified by the 
fact that the system or measurement noise is typically 
caused by a number of sources. It can be shown 
mathematically that when a number of independent random 
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variables are added the resultant effect is very nearly a 
Gaussian probability density even though the individual 
densities are not Gaussian [Ref. 9l* 

Under the above mentioned assumptions of a purely 



white Gaussian noise, the first two moments specify the 
entire shape of the density describing the noise, and the 
mathematics of the problem are greatly simplified. 
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II. KALMAN FILTER EQUATIONS 



A. GENERAL 

The design of a Kalman filter and especially the 
integrated I.N.S. Kalman filter design requires extensive 
computer simulation. This chapter of the work is a 
presentation of the equations which are required not only 
for the mechanization of the filter but also those which are 
necessary to simulate the dynamics of any user (aircraft, 
missile) . 

The principal tool used for the solution of this 
specific and other similar problems is the very common 
method of covariance analysis. It is known that the 
covariance is a measure of the uncertainty in the knowledge 
of the true values of the state vector components. In this 
work, as the covariance matrix is concerned, the off 
diagonal terms are assumed to be zero initially and initial 
conditions on the diagonal elements are arbitrarily taken. 
The covariance matrix of both the system and the filter are 
propagated forward in time by numerical integration 
techniques . 

The adjustment of the values of the state variables, to 
those of the best estimate obtained with the Kalman filter , 
is achieved when a control is applied to the system after 
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the specified update time is reached and the best estimates 
of the states have been determined. The square root of the 
individual diagonal elements of the system covariance matrix 
(RMS values) are plotted as a function of time to provide 
the performance of the filter. For this study the plots are 
also utilized to determine the error contribution associated 
with each modeled error source. Furthermore, the error 
statistics are propagated which means that the standard 
deviation of the noise value is supplied whenever a noise is 
required . 

One attribute of covariance analysis is, that under the 
assumptions stated in Chapter One, i.e., the linearity and 
white Gaussian noise, the covariance is independent of the 
actual measurement values and can be computed through 
generating a sample sequence of measurements. And as a 
matter of fact this method is easier to handle and work with 
than the corresponding Monte Carlo type simulation. 



B. SYSTEM MODEL EQUATIONS 

The differential equations that describe how the 
inertial nav iga tor errors propagate with t ime are the basic 
equations used in this process. These equations are 
formulated in a set of first order, linear differential 
equations, driven by white Gaussian noise for the reasons 
described previously in this work. 
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Linear measurements corrupted also by white Gaussian 
noise are made upon the actual system variables. It is 
furthermore assumed that the equations which represent a 
detailed model of the system are of the form: 



— s s— s 3— s 



( 1 ) 
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w, is an m vector of white noise inputs with the 



E[w( i) w( j)^] 



Qg( i) for i = j 



( 2 ) 



0 for i j 

where the indices i and j are instants in time. 

The observations which are obtained from external 
references and in our case of study from the G.P.S. can be 
described by the following linear measurement vector 
equation : 



Zc = H^x^ + v„ 
— s s— s — s 



( 3 ) 
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where : 



z is a q vector of measurements 



is a q X matrix of measurements 

Vg is a q vector of white noise inputs with the 
characteristics of zero mean and variance: 



E[v(i) v( j)^] 



R^(i) for i = j 

0 for i j 



( 4 ) 



A further assumption for the study is that the 
w and the measurement noise v are uncorrelated 
i . e . , 



system noise 
for all time, 



E[w(i) v( j)^]* = 0 



for all i,j 



( 5 ) 



C. FILTER EQUATIONS 



The equations discussed above are assumed to be a 
complete and accurate mathematical description of the G.P.S. 
aided inertial navigation system dynamics and measurement 
equations for the purpose of simulation. They also 
constitute a set of equations which would be utilized in the 
design of a fully optimal Kalman filter. 

In our case of study as also in general a suboptimal or 
reduced order filter design is obtained by reducing the 
dimension of the state vector due to the computational 
burden of the fully optimal filter. The states that are 
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eliminated are those that affect the accuracy the least of 
the mathematical description of the aided-I.N.S. The 
designed suboptimal filter can be implemented with the on- 
board computer (aircraft or missile). 

The equations which describe the suboptimal filter are 
of the form: 

• 

iif ~ ‘^f— f *^f— f 

where 

is an ri2 vector 

is an CI2 ^ 1^2 filter dynamics 

Gj-. is an x m2 matrix of gains 

w^ is an m vector of white noise inputs with the 
characteristics of zero mean and variance: 

^ Q^(i) for i = j 

E[w^(i)Wj:.(j) ] = (6) 

0 for i j 

The equation for the filter measurement is : 

—f - ^fiif If 

where : 

is a q vector 
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is a q X q matrix of measurements 

is a q vector of white noise inputs with the 
characteristics of zero mean and variance: 

rr, Rf.(i) for i = j 

ECv^(i) v^( j) M = (3) 

0 for i j 

The filter propagation and update equations based on tne 
above models are then given below. 

At measurement times (update): 

(9) 

P^^ = P^“ = (10) 

and between measurements (extrapolate): 

If = fflf ^^2) 

P^ = F^Pf + (13) 

where : 

is an r \2 vector denoting the best estimate 
P^ is the covariance matrix of the filter 
is the matrix of Kalman gains 
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is a q vector of the actual values of the 
measurements taken 



+ superscript indicates the time instant just after 
update 



- superscript indicates the time instant just prior 
to update 



T superscript denotes the 
super scr ipted . 



transpose matrix or vector 



The filter subtracts from the actual taken measurement 
z the best prediction of its value before the actual 
measurement is taken, i.e., the value of This 

X 

difference is then passed through an optimal weighting 
matrix and used to correct " > the best prediction of 
the state at the time instant before the measurement is 
taken. This process gives the best estimate after update. 
This estimate is propagated to the time of the next 
measurement sample according to equations (12) and ( 13 )- 
The above recursive relationships are solved based on 
initial conditions of an assumed Gaussian density which 
describe the a-priori knowledge of the state as: 



x(0) 







and P(0) = Pq 



( 14 ) 
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The Kalman filter conditioned on the actual measurements 



taken, propagates the conditioned probability density of the 
desired states. The probability density function of a 
Gaussian noise amplitude takes on the shape of a normal 
bell-shaped curve. The assumption of Gaussian noise 
amplitude is well justified by the fact that a system or 
measurement noise is typically caused by a number of small 
sources and according to the Central Limit Theorem it can be 
shown mathematically that when a number of random variables 
are added together, the summation is a random variable whose 
density is nearly a Gaussian probability density, regardless 
of the shape of the densities of the individual random 
variables. Furthermore, the use of Gaussian densities makes 
the mathematics easier to handle and tractable. It is known 
that a Gaussian density is completely determined by its 
first and second order statistics, i.e., the mean and the 
variance. Thus, the Kalman filter, which propagates the 
first and second order statistics, includes all information 
contained in the conditional probability density mentioned 
above [Ref. 2: pp . 3-9]* 

The mean of a density function or its expectation u, is 
defined as 




xf ( X ) dx 



(16) 



— . 00 
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and it is interpreted as the weighted average of the values 
of X, using the probability density function f(x) as the 
weighted function. All the Gaussian white noise inputs in 
this study are assumed to have zero mean. 

The variance of a density function^or the square of the 
standard deviation o, is defined as: 



Var[x] = a~ 



J ( X - u)"f ( x) dx 



(17) 



and it is interpreted as the weighted average of the values 
of (x - w) ; thus, a“ is a measure of the density spread and 
a direct measure of the uncertainty since the larger a is, 
the broader the probability peak is, spreading the proba- 
bility weight over a larger range of x values. For the 
example of Gaussian density, 68.3"^ of the probability weight 
is contained within the band of a units to each side of the 
mean u, which represents the area under the normal bell- 
shaped curve between the values of -a and +a and 95.4% of 
the probability weight is contained between the values of 
-2 0 and +2a. Since in this study vectors are used instead 
of scalars, the above equations (16) and (17) which give the 
first and second order statistics respectively for the 
scalar case, must be extended for the vector case as 
follows : 
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£ f (jc) dx ^ • 



.dx^ 



(18) 



n 



Cov[jc] = P = / ... / (x - u)(x - u)^f( x)dx^ . . .dx,^ (19) 

^ 00 * CO 
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III. I.N.S. ERROR ANALYSIS 



A. GENERAL 

One general approach to determine the navigation error 
caused by individual sources of error is to simulate the 
inertial-navigation-systam nonlinear equations and sources 
of error and compare the navigation outputs witn tne 
simulated true position — the difference being the navigation 
error. However, this is not the approach used here. 

Assuming that tne position errors are small compared with 
earth radius, that the velocity errors are small compared 
with orbital velocity, and that the alignment errors are 
small compared with 1 radian (or 3^37 arc-min) it can be 
demonstrated that the propagation of errors in an inertial 
navigation system is very accurately governed by a set of 
corresponding linear differential equations. Therefore, 
most inertial-navigation-system error analyses are conducted 
working directly with a set of linear error differential 
equations . 

Sets of error equations have been developed for various 
I.N.S. configurations in the context of a particular 
application. As a consequence, many sets of I.N.S. error 
equations have been developed for the various broad classes 
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I 












of mechanization such as local level, space stabilized, 
wander azimuth, free azimuth, strapdown , etc. 

The choice of navigation error variables for analysis 
has often followed from the coordinate system used in the 
navigation equations or that implied by the physical I.N.S. 
platform orientation. Regardless of the differences in the 
sets of equations the fact is that I.N.S. error propagation 
is to a large extent, completely independent of system 
tuechani zation . Britting [Ref. 2] has shown that the basic 
error differential equations for any I.N.S. may be written 
in standard coordinates, regardless of the physical 
mechanization or internal navigation variables. Further- 
more , the unforced (homogeneous) portion of these 
differential equations is, under certain very broad 
assumptions, identical for any arbitrarily configured 
terrestrial I.N.S. 

The error equations presented in this chapter follow the 
philosophy of [Ref. 2] including the choice of north-slaved 
coordinates for the error variables and the identification 
of the unforced (homogeneous) portion of the differential 
equations that is independent of system mechanization. Two 
of the differences that are noticeable are: 

1) The error equations are written as a system of nine 
first order equations (and further on reduced to seven) 
rather than three second order plus three first order 
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equations. The first-order form is the state space repre- 
sentation of the error equations used in modern estimation 
theory. 

2) The form of altitude compensation assumed in 
[Ref. 2] is not found in most inertial navigators. Gravity 
is assumed computed as a function of the inertial system 
indicated position. 

B. GENERAL I.N.S. ASSUMPTIONS 

The assumptions pertaining to the general error 
differential equations are broad enough to encompass all of 
the important I.N.S. configurations [Ref. 2]. 

1) Three accelerometers are available to measure the 
specific force vector. The equations for a two-accelero- 
meter local level system are the same provided the inertial- 
altitude and vercical-velocity equations are deleted. 

2) The accelerometers are mounted on a platform whose 
angular velocity is either controlled (as with gyro- 
stabilized gimbaled platforms) or is measured (as with the 
strapdown systems) . 

3) The system's indicated velocity vector and three- 
dimensional indicated position are obtained by integration 
of the gravity field compensated specific force measure- 
ments, using a correct set of differential equations. 

4) A model of the earth's gravity field is used to 
calculate the gravity vector as a function of the system- 
indicated position . 



42 



5) External altitude information and other navigation 
measurements may be used to update the inertial-navigation- 
system indicated position, velocity, and attitude. 

6) A computer is availabe to process the navigation 
information and the computation errors are either negligible 
or may be treated as equivalent instrument uncertainties. 

7) Both the mechanical coordinate frame (the frame 
tracked by the platform) and the computacion frame (the 
frame to which the specific force measurements are 
transformed for velocity and position integration) are 
arbitrar y . 

C. LOCAL-LEVEL TERRESTRIAL NAVIGATOR 

Many inertial navigators do use a local-level coordinate 
system for the velocity and position integration. The 
local-level terrestrial navigation system physically instru- 
ments the local geographic coordinate frame. The platform 
axes are commanded into alignment with the local north-east- 
down coordinate system. 

The local-level terrestrial system is undoubtedly the 
most successfuly of all inertial-navigation-system configur- 
ations. The class of local-level systems today constitutes 
the majority of operational inertial navigations systems. 
Since this system is described in detail in [Ref. 2] here 
are listed some of its advantages: 
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1) The computation of gravity components is greatly 
simplified. In fact some navigators use zero for both 
horizontal components of gravity. 

2) Some inertial navigators have no vertical accelero- 
meter (this is the case in the present study) and do not 
mechanize a vertical channel. The horizontal velocity and 
position equations of a local-level set are appropriate for 
such a nav ig ator . 

3) The well-known altitude and vertical velocity 
instability of a pure inertial navigator must be stabilized 
by means of an external altitude reference. But a local- 
level set of variables includes altitude and vertical 
velocity explicitly. The altitude stabilization equations 
therefore can be simplified. 

4) The calculations required to provide navigation 
outputs and displays in geographic coordinates are 
simplified . 

D. THE TOO ACCELEROMETER LOCAL-LEVEL SYSTEM 

Many inertial navigators have only two accelerometers 
after the vertical accelerometer has been eliminated. The 
system is composed of a three axis inertial platform, two 
accelerometers which are nominally orthogonally mounted in 
the instrumented east and north directions and a computer 
which performs the necessary navigational computations 
[Ref. 2]. 
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The north and east gyros are respectively connected with 
the instrumented north and east accelerometers at the signal 
level and the gyros are torqued at a rate proportional to 
the vehicle's longitude and latitude rates so that the 
platform can maintain its axes aligned with geographic axes 
since the vehicle carrying the navigation system is assumed 
to move freely over and above the earth. The accelerometer 
outputs provide these required torquing signals which must 
be so compensated that gyro command can be obtained as a 
function of only the north and east velocity rates. 

Such a two-accelerometer local-level I.N.S. has seven 
state variables: two of position, two of velocity, and 

three of platform alignment. 

1 . Error Model Equations 

The general model of local-level inertial navigation 
systems is given by the following matrix equations: 

AX = q (20) 



where 

A = system characteristic matrix, same for all I.N.S. 
configurations 



q = forcing vector of inertial system errors 

T 

or ^ = [ Q 1 > Q2’ ^3’ ' ^ 6 ’^ 



( 21 ) 
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X 



error state vector of attitude, position, and 



velocity errors 



or ^ sL, 6l, 5ii] 



T 



( 22 ) 



It is quite important to emphasize that a computer simula- 
tion program developed in accordance with the above equation 
(20) is valid for all possible I.N.S. configurations, space 
stabilized, strapdown , wander-azimuth and not only for the 
local-level one. Both the coordinate frame mechanized by 
the inertial instruments and the computation frame are 
completely arbitrary. It is only the forcing function, q, 
which depends on the system configuration through the angu- 
lar velocity and orientation of the inertial instruments. 



vector-differential equation the error state of the I.N.S. 
is defined as [Ref. 2]; 



and for the case of this study for two-accelerometer local- 
level I.N.S. system this reduces to: 



where the seven basic I.N.S. errors are: 

^N’ " north, east, down platform tilt errors 



In order to rewrite this equation as a first order 




(23) 




(24) 



6L , 61 



latitude, longitude position errors 



46 



1 



^57 


= 1 


^62 


= -fo/^ 


^63 


= fg/r 


^64 


= -1(1 + 2w^^)cos 2L 


^67 


= -X sin 2L 


^71 


= fp/r cos L 


^73 


r -f[^/r cos L 


^74 


= 1 tan L + 2XL 


■^76 


= 2x tan L 


F 

77 


= 2 L, tan L 


e = 7 X 


1 error vector matrix 



where 



£ “ 


Sg, Cq, 6L, 6l, 6L> 51] 



= 7 X 5 forcing matrix, with non-zero elements 





II 


^22 


= 1 


^33 


II 



47 



6L, 51 = latitude, longitude rate errors 

The above statements allow equation (20) to be written as: 

€ = + Ge^ (25) 

where : 
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7 X 7 error matrix with non-zero elements 
-X sin L 
L 

-X sin L 
cos L 
X sin L 
X cos L 
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-X cos L 
- sin L 
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1/r cos L 



and q = 5 X 1 forcing vector matrix 



where 



q 




q p > 3 ’ *4 ^ > 




(23) 



and neglecting both gyro and accelerometer non-orthogonalioy 
errors the forcing functions are comprised of 10 I.N.S. 
component errors as below: 
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where for our study of two accelerometer local-level system 
the platform-to-navigation transformation matrix and the 

I.N.S. platform torquing rate below: 
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Finally the forcing vector matrix is modeled for our case 
as : 



" 1 




— "■ 
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2 . Error Equations Solution 

The solution of the differential equations 
represented by equation (25) gives the error response for 
the two-accelerometer local-level navigator for arbitrary 
vehicle motion within the constraints implied by a "first- 
order" analysis. Since the coefficients of the differential 
equations are time varying the analytic solution of equation 
(25) would be quite tedious and require the user of a 
computer program to generate flight profiles. In our study 
specific cases are examined so that the coefficients of the 
differential equations can easily be calculated. 

The specific cases which are examined are the 
following three: 

•• «• «• 

• • • • 

1) Stationary case, where \ = and L=L = l = l = h=h = 0 

2) Easterly flight at 600 "ft/sec" or 355-5 "knots," 
where x ~ 1.557 x 
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3) Westerly flight at 600 "ft/sec'» or 355.5 "knots,*' 

where X = 0. 442 x u . 

le 

Writing the equation (25) in terras of the error 
states explicitly we have the following system of differen- 
tial equations which must be solved simultaneously: 

x(1) = -x(sinL)x(2) - X(sinL)x(4) + (cosL)x(7) + q -] 

x( 2 ) = x(sinL)x( 1 ) + x(cosL)x(3) - x(5) + q 2 

x(3) = -X(cosL)x(2) - X(cosL)x(4) - (sinL)x(7) + 4^ 

x(4) = x(6) (32) 

x(5) = X (7 ) 

x( 6 ) = p x( 2 ) - x(sin2L)x(7) + 14 

= - r(Jsr y >“'> ^ 2;(tanL)x(6) . qj 

The values of the parameters used for the computer 
simulation are given in Table II. 

a. Constant Gyro Drift Errors 

Letting constant gyro drift be the sole error 
source in the I.N.S. system where the constant gyro drift 
rates (u)uj^, (u)o)£, (u)a)j^ are associated with the north, 
east, and azimuth gyros respectively, computer simulation 
verfied the following. 
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TABLE II 



COMPUTER SIMULATION PARAMETERS 



Geographic Latitude, L 


45°L 


Terrestrial Longitude, 1 


0*^ (Greenwich) 


Earth Rate , 

Constant Gyro Drift 


0.2618 Rad/hour 


• North, (u)ujj^ 


1 meru ( 1 ) 


• East , ( u ) 


1 meru 


• Azimuth , ( u) to^ 
Constant Accelerometer Bias 


1 meru 


• North, (u)fj^ 


200ug (2) 


• East , ( u) f c 
Initial Platform Tilt 


200ug 


• North, £j^(0) 


0. 14 mrad 


• East, e-(0) 


0.14 mrad 


• Azimuth , £q( 0) 


0.14 mrad 


Initial Latitude Error, 5L(0) 


0.17 mrad (3) 


Initial Longitude Error , 51(0) 


0.17 mrad 


• 

Initial Latitude Rate Error, 5L(0) 


0.34 mrad (4) 


Initial Longitude Rate Error, 51(0) 


0 . 34 mrad 


Stationary Flight 


X = 


Constant Easterly Flight at 600 ft/sec 


X = 1 . 557xu^g 


Constant Westerly Flight at 600 f t/sec 


X = 0.442X(u^g 


Gravitational Acceleration Constant 


g = 32.2 ft/sec 


Earth Radius 


r = 6,37 8.3 :<m 


1. The RMS gyro drift of 1 meru corresponds 


to 0 . 015 °/hr or 



2 . 

3. 

4. 



0.2618x10“-’ rad/hr. 

The RMS accelerotnete^ bias of 200ug corresponds to 
0. 4356x10“^ rad/(hr)‘^. 

The RMS position error of 0.17 mrad corresponds to 1085tn 
or 0.586 arc min . 

The RMS velocity error of 0-3^ mrad corresponds to 
2 ft/ sec . 
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For the stationary case we observe in Figures 4 
through 10 that for the north and east level errors, ej,^ and 
Eg, the Foucault modulation is an effect of first-order in 
contrast with the latitude, longitude and azimuth errors 5L, 
5l, and eq respectively where the Foucault modulation has 
only a second-order effect [Ref. 1]. These computer 
solution results suggest that it would be convenient, for 
design purposes, to neglect the Foucault modulation since 
the equations we obtain then are easily solved and give 
solutions with approximately the same amplitude information 
for latitude and longitude which are of primary importance 
for navigational purposes while the relatively poor 
information of one level errors is of secondary importance. 

We further note from these computer graph 
solutions that the effect of the Foucault terms in ohe error 
equations system, is to modulate the Schuler oscillations at 
a frequency given by the local vertical projection of earth 
rate, namely x sin L, which corresponds to a period of 

33-9 hours for the selected latitude L = 45°* 

For the case of constant easterly flight at 600 
ft/sec the results are given in Figures 11 to 17. Compari- 
son with the curves for the stationary case indicates that 
the lowest modulation frequency has increased from i 
to X z 1.557 X and the space rate period is 10.9 hours 

while the Foucault modulation now occurs with a period of 
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about 21.8 hours instead of the 33*9 hours period for the 
stationary case . 

Another important feature revealed by the above 
comparison is that the azimuth and latii^ude errors are 
reduced from the corresponding for the stationary case by a 
factor of 1.557 which represents the ratio for this 

case . 

For the responses to the north and azimuth zero 
drift, (u)ujq and (u)uj£, the vehicle motion appears to have 
little effect on the error growth in the cases that exhibit 
a longitude error which grows with time. 

The level errors in response to level gyro drift 
are seen to remain unchanged while the level error response 
to azimuth gyro drift, (u)uq) is seen to emerge from the 
computer noise having a peak value of 2.3 rad/meru (Figure 
13) . 

The longitude error in response to azimuth gyro 
drift which was bounded for the stationary case is now 
reduced by the factor or by 1.55 while the latitude 

and longitude rate error magnitudes are unaffected by the 
vehicle motion . 

Finally for the westerly flight case it is 
verified that the level errors remain unchanged without the 
effect of Foucault modulation, but the latitude, longitude, 
and azimuth errors grow approximately in proportion to the 
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time-drift rate prodct. The computer solution graphs for 
the westerly flight case are shown in Figures 19 through 24. 

Similar results are found for the cases of east 
and azimuth gyro drift but they are not included here due to 
large amount of graphs. 

b. Accelerometer Bias Errors 

Considering the accelerometer bias as the sole 
error source computer simulation of equations (25) shows the 
following result of the effects of the north and east 
accelerometer bias, (u)fj^ and (u)f£ respectively, on the 
navigation and level errors. 

For the stationary case the results are shown in 
Figures 25 through 31* note that the Schuler mode is 

predominant since the accelerometer bias directly excites 
the relatively high gain level loops and that the Schuler 
oscillations are modulated at the Foucault mode frequency of 
33*9 hours per cycle. The maximum values for the navigation 

_7 

errors proved to be in the range of 20 x 10 rad/200 yg for 

the latitude error and 1.4 x 10“ rad/200 ug for the 

longitude error. 

We notice as for the constant gyro drift case 
that we can neglect the effects of Foucault modulations as 
first-order ones and proceed in the solution of the 
resulting equations more easily obtaining almost the same 
approximate amplitude information of the navigation and 
level errors . 
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For the case of constant easterly flight at 600 
ft/sec the computer simulation graphs are shown in Figures 
32 through 38 . We observe here that again the Foucault 
modulating frequency has increased by a factor of 1.55 which 
corresponds to the ratio Nevertheless we see that 

the error sensitivities remain unchanged with the previously 
explained stationary case. 

Figures 39-45 snow computer solutions of the 
navigation and level errors for the case of westerly flight 
at 600 f t/sec . 

We see now that the Foucault modulating 
frequency has decreased by a factor 0.442 correspond ing to 
the ratio this case and once again we can proceed 

to the solution of the equations without considering the 
Foucault terras, especially for design purposes. kn easy 
extension of the above observations is that for the limiting 
case when the terrestrial longtitude rate, 1 , in a western 
direction is equal to the earth rate, then the Foucault 

modulation disappears completely leaving a pure Schuler 
oscillation . 

c. Initial Condition Errors 

The results on the navigation and level errors 
due to effects of the initial conditions are now presented 
accompanied by only the most important graphs of the 
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computer simulation since the total amount of graphs is too 
large to be included in this study. 

For the stationary case selected graphs 
representing the level errors for an initial north, east and 
azimuth level error of 0.14 rarad or 0.438 arc-min are shown 
in Figures 46-47, 43-49, and 50-51 respectively. 

Figures 52 through 57 present the level and 
navigation errors for an initial latitude rate error of 2 
ft/sec corresponding to 0.34 mrad/hour while figures 53 
through 63 show the associated errors with an initial 
longitude rate error of the same amount. 

There is no need to include any graphs for the 
resulting errors due to initial longitude error since by 
inspection of the error differential equations we can see 
that longitude is uncoupled from the other computation loops 
so that an initial longitude error holds constant and no 
other error becomes non-zero . 

We discussed up to now the errors of a pure 
I.N.S. system. In the next chapter we proceed with the 
consideration of the combined I . N . S. /G . P . S. system and the 
results we achieved after computer simulation. 
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IV. I .N .S ./G .P .S . SYSTEM MODEL AND EQUATIONS SOLUTION 



In order to apply the Kalman filter equations discussed 
in Chapter II a reference system model which is a good 
approximation to the real world dynamics is needed. 

In this chapter we outline the reference I. N. S./G . P. S. 
system equations selected for this study. First we are 
defining the error states incorporated in the system model 
along with their assumed initial conditions values. Next we 
discuss the modeling of the I.N.S. plant error states and 
finally we present the equations for the integrated 
I . N. S./G . P. S. system together with their simulated computer 
results . 

A. SYSTEM MODEL 

1 . State Variable Definition and Initial Conditions 
In Table III we present a listing of the state 
variables utilized in the reference system model. The 
initial conditions on the I.N.S. error states are highly 
arbitrary and the selected values are similar to those used 
in other unclassified studes [Refs. 4 , 51 . 

For the initial conditions on the gyro error states 
and the accelerometer, the values are selected for a typical 
inertial navigation system of one nautical mile per hour 
class . 
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After the above definition of the initial conditions 
we have by the same time specified in a complete w.ay the 
initial covariance matrix P(0), since its diagonal elements 
are the squared values of the given RMS initial conditions. 
The remaining off-diagonal elements of the initial covar- 
iance matrix are assumed to be zero initially. 

Furthermore the propagation of the linear variance 
equation (13) requires an additional knowledge of the two 
matrices F and Q* where; 

Q* = GQG^ (33) 



where 

G is the forcing input matrix 

Q is the input noise covariance matrix 

and the F matrix is the same as in Equation (25) and which 
has been used in the previous chapter for the inertial 
navigation system error equations solution and computer 
simulation . 

For the Q* matrix the only non-zero elements are all 
diagonal and we will denote these from now on as where 
the subscript i denotes the row and column of the value. 

For example, indicates that this is the value which 
belongs to the intersection of the 3rd row and the 3rd 
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column in the Q* matrix and cor 
input on state variable number 
These non-zero elements 
matrix are five, corresponding 
and 7 according to the notation 
2. Plant Error States 
The fo 1 1 o wi 
level errors, X and 
errors, constitute 
equations of these 
dynamic response of 
system . 

Tnere are v 
these error states, 
will use again the 
matrix F given in e 
local- lev el two-acc 
con f ig ur ation . 

B. EQUATIONS SOLUT 
Using the defin 
write the following 

_x( t) = F jc( 

z( t) = H x( 



responds to a white noise 

3. 

in the reference system Q* 
to state numbers 1, 2, 3> 6 
of Table III. 



ng seven states, North, East, and Azimuth 
Y position errors, X and Y velocity 
the plant error states. The differential 
states describe the natural unforced 
the errors in the inertial navigation 

arious models for the implementation of 
As we did in the previous chapter we 
Pinson error model described by the 
quation (25) for our specific case of the 
elerometer inertial navigation system 



ION AND COMPUTER SIMULATION 
itions described in previous pages we can 
equations for the error states: 



t) + G w( t) 
t) + v(t) 



(34) 
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TABLE III 



I.N.S./G.P.S. SYSTEM STATE VECTOR DEFINITION 

Error 



State 


Syrabol 


Definition RMS 


1 Initial 


Cond i tion 






I.N.S. PLANT ERROR STATES 




1 . 




North atttiude error 


0. 14x10"^ 


Rad 


2 . 


£ 

il 


East attitude error 


0. 14x10"^ 


Rad 


3. 


^ D 


Azimuth attitude error 


0.14x10'^ 


Rad 




u 






(1) 


4. 


5L 


Y position error 


o 

X 

o 

1 

L 


Rad 


5. 


51 


X position error 


0. 17x10"^ 


Rad 


6. 


5L 


Y velocity error 


0.34x10"^ 


(2) 

Rad/hour 


n 

{ • 


51 


X velocity error 


0. 34x10"^ 


Rad/hour 






I.N.S. ERROR SOURCES 






8. 


( U ) u) 


North Gyro Drift 


1 meru 


(3) 


9. 


(U)ojg 


East Gyro Drift 


1 meru 




10. 


( w ^ 


Azimuth Gyro Drift 


1 meru 




1 1 . 


(u)f, 


North Accelerometer Bias 


200xl0"^g 


(4) 


12. 


(u)fE 


East Accelerometer Bias 


200x10’°g 





1. The RMS position error of 0.17 milliradians correponds to 
1085m or 0.586 arc rain. 

2. The RMS velocity error of 0.3^ raillirad/hour corresponds 
to 2 ft/ sec . 

3. The RMS gyro drift of 1 raeru corresponds to 0.015 /hr or 

261.8x10"° rad/hour. . 

4. The RMS acceleroraetec bias of 200x10"°g corresponds to 
0.4356x10"^ rad/Cnr)^. 



62 



where 



£(t) =7x1 error state vector 

- ^E’ ^1] 

£ =7x7 Pinson error model matrix 

as described in equation (25) 

G =7x5 input forcing matrix as described 
in equation (25) 



_z(t) =2x1 vector of measured states 
H =2x7 measurement matrix where 



H 



0 0 0 1 0 0 0 
0 0 0 0 1 0 0 



denoting that we have available measured information for ohe 
X and Y position error . 



w(t) = 5 X 1 forcing vector assumed to be white 
Gaussian noise 



and 



_v(t) =2x1 vector of measurements noise assumed to 
be white Gaussian. 



Using the feedback configuration of the Kalman filter we 
can write the following equations: 
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^(t) = F £(t) + G w(t) + X[£ - H £(t)] 

and 

P(t) =F£+PF^+GQG^-PH^R“^HP 
where 

^(t) = the covariance matrix 

K = Z. the Kalman filter gains matrix 

R = the measurement noise covariance matrix 

In order to achieve numerical results via computer 
simulation we write the predicted error states of our system 
in explicit form as below with the help of Table IV in which 
the system states and their corresponding symbols are 
defined : 

x( 1 ) = -X(sinL)x(2) - X(sinL)x(M) + (cosL)x(7) + AA + 

+ K 3 _^[x( 8 ) - x(4)] -h K^ 2 tx( 9 ) - x(5)] (35-a) 

x(2) = x(sinL)x( 1 ) + x(cosL)x(3) - x(6) + 3B + 

+ K 2 ^[x( 8 ) - x(4)] + K 22 [x( 9 ) - x(5)] (35 -d) 

x(3) = -X(cosL)x(2) - X(cosL)x(4) - (sinL)x(7) + S3 + 

+ K^^[x( 8 ) - x(4)] + K^ 2 t^( 9 ) - x(5)] (35-0 

^4) = x( 6 ) + Ki^^[x( 8 ) - x(4)] + K42tx(9) - x(5)] (35-d) 



(35) 



(36) 
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TABLE IV 



COMPUTER SIMULATION VARIABLES AND CONSTANTS 



x( 1 ) 


= 


^ N 




= 


North attitude error 


x(2) 


r 


^E 




= 


East attitude error 


x(3) 


= 


"d 




z 


Azimuth attitude error 


x(4) 


= 


5L 




z 


Y position error 


x(5) 


= 


51 




z 


X position error 


x(6) 


z 


6L 




z 


Y velocity error 


x(7) 


= 


61 




= 


X velocity error 


x(8) 


z 


O 




= 


G.P.S. Y position error measurement 


x(9) 


= 


O 




= 


G.P.S. X position error measurement 


X (10) 


= 


“-t 




= 


True Y position error 


x( 1 1) 


= 


“t 




z 


True X position error 


x(12) 


= 


W 






Input white Gaussian noise 


x(13) 


z 


V 




= 


Measurement white Gaussian noise 


A 


Z 


1 — 1 
c 
e 




z 


North Gyro Drift Variance 


B 


= 


[(U)0J 




z 


East Gyro Drift Variance 


S 


z 


[(u)o) 




Z 


Azimuth Gyro Drift Variance 


D 


= 


C(u)f 


n’' 


= 


North Accelerometer Bias Variance 


E 


z 


C(u)f 




z 


East Accelerometer Bias Variance 


F 


z 


«11 




z 


Measured Y position error variance 


G 


= 


^22 




z 


Measured X position error variance 


AA 


= 


x(12) 


•A 


z 


White noise like North Gyro Drift strength 


BB 




x( 12) 


•B 


- 


White noise like East Gyro Drift strength 
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r 




TABLE IV (CONTINUED) 



S3 




x( 12) *3 




White noise 
strength 


like 


Azimuth Gyro Drift 


DD 


- 


x( 12) *D 


- 


White noise 


like 


North Accelerometer 










bias strength 




EE 


- 


x( 12) *E 


- 


White noise 


Like 


East Accelerometer 










bias strength 




FF 


= 


x( 13) *F 


= 


White noise 


like 


Y-position error strength 


GG 




x( 13) *G 


= 


White noise 


like 


X-position error strength 


k 




g/r 




constant 
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x(5) = x(7) + K^^CxCS) - x(4)] + K^2^x{9) - x(5)] (35-e) 

/V 

x(6) = kx(2) - x(sin2L)x(7) + DD + 

+ Kg^Cx(8) - x(4)] + - x(5)] (35-f) 

x(7) = - x(1) + 2x(tanL)x(6) + EE + 

+ K^^[x(8) - x(4)] + - ^^5)] (35-g) 

Assuming the input forcing vector as white Gaussian 
noise whose strength is related to the value of the variance 
of each input error source (the corresponding one) computer 
simulation was proceeded in the following way. 

First with the help of the RICATI FILTER computer 
program available at the NFS W.R. Church Computer Center we 
solved the corresponding for our study Ricati equation of 
covariance propagation obtaining the Kalman filter gain 
matrix. The data we used to run the above program are 
provided in Tables III and V. A listing of the data 
formulation for the RICATI FILTER program is given in 
Appendix D. 

The calculated from the above program values of the 
Kalman filter gains for a processing period of four hours 
are given in the following Table VI. 

Additional runs of the above program have been contacted 
for processing periods up to 36 hours and it has been 
observed that the Kalman filter gains reach a steady state 
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TABLE V 



NUMERICAL VALUES FOR RICATI PROGRAM ^ 
F Matrix (7 x 7 ) 



12 


= 


-0.1851 


Fi 4 = 


-0. 185 1 


F17 = 


0.707 


21 




0.1851 


"23 = 


0.1851 


^26 - 


-1.0 


32 




-0. 1851 


F34 = 


-0. 1851 


F37 = 


-0.707 


46 


= 


1.0 










57 


= 


1 . 0 
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r 


19.92 


^67 ^ 


-0. 2613 






71 


r 


-28.175 


^76 = 


0.5236 







All other elements are zero. 



Matrix (5 x 7 ) 

" ^22 ~ ^33 ' 

G 45 = 0.0915 

= 0.1294 

All other elements are zero. 

H Matrix (2 x 7) 

^14 ~ ^25 ~ 

All other elements are zero. 



For the used values of elements in the matrices the 
results will be given having units the appropriate for each 
case combination of radians and hours. 
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TABLE V (CONTINUED) 



Q Matrix (5 x 5) 

Qn “ ^22 ~ *^33 ■ to. 0002618 rad/hr]^ = 0.7x10 
= [0.0004358 rad/(hr)^]^ = 0.19x10"^ 

All other elements are zero . 



R Matrix (2 x 2) 

^11 = ^22 ■ ‘3*52x10"'^ [rad]^ 
^12 ■ ^21 “ 



P(0) Matrix (7 x 7) 



P^^(O) = Pp2t3) = ?^^(0) = [0.00014 rad]- = 0. 
P^^(O) = ^55(0) = [0.00017 rad/hr]- - 0.3x10“^ 
P.g(O) = P^^(O) = [0.00034 rad/(hr)^]^ = 

0.115 X 10”° [rad/(hr)^]^ 

All other off-diagonal elements are zero. 



”^[rad/hr]^ 

[rad/(hr)^]^ 



1 x1 0” ^ [ rad ] ^ 
[rad/hr ] - 
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TABLE VI 



KALMAN FILTER CASE 

K^^ = 0.034756672 

K2-, = -1.36921265 
K3. = 1.56555360 

Kj^^ = 3.4402061 1 

= 0.0907281520 

= 5.92162434 

= 1.04245351 



FOR A 4 HOUR PROCESS 



Ki2 = 


1 . 17201 123 


II 

C\J 

CM 


-0.0401669175 


K32 = 


-2.86193609 


II 

CM 


O.O90728I5O 


II 

CM 


4.31037639 


^62 = 


-0.339257459 


II 

CM 


9. 2937881 1 
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condition for which their values are not very much different 
from those achieved for a 4 hour processing period. 

So in the following calculations and computer simulation 
we have used the values of the Kalman filter gains resulted 
for a 4 hour simulation period. 

Having available the values of Kalman filter gain matrix 
elements which are used to multiply the residuals in the 
appropriate equations in order to achieve the predicted 
error states of the integrated I . N . S . /G . P . S . system, the 
appropriate program has been formulated in order to solve 
the error differential equations ( 32 ) described above, with 
the use of the available routine INTEG2S slightly modified 
for accurate evaluation and plot of the error state 
variables . 

The simulation results for the I . N . S. /G . P . S . system 
operation for a period of four (4) hours are presented below 
in Figures &4 through 93 . We can easily observe that all 
the state error variables of the combined I. N.S./G. P.S. 
system are damped out and the resulting value of the errors 
after a period of one (1) hour is small enough so that the 
operation of our system model can be characterized as 
successful . 

Specifically in Figure 64 we observe that the north 
level error of our system is dropped down to 0.025 milli- 
radians after one (1) hour even if the starting initial 
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condition value was 0.14 milliradians . Furthermore at the 
end of a 4 hour period the error has been diminished to the 
value of 0.005 milliradians which is subject to error 
reduction by a factor of 28. 

In Figures 65 and 66 the East and Azimuth attitude 
errors are presented respectively where similar as with the 
north attitude error observations occur. 

The Y position error behavior is presented in Figure 67* 
There we can see that even if we started from an initial 
condition error of 0.17 milliradians (or 3256 ft) after one 
hour processing the error has been diminished to only 0.028 
milliradians (or 536 ft) and furthermore after a four hour 
period this error drops down to 0.01 milliradians (or 
191 .5 ft) . 

The X position error damping out seems to be more 
attractive since from Figure 68 we can see that after one 
hour the error drops down to 0.020 milliradians (or 383 ft) 
and at the end of a four hour period the error is diminished 
to 0.0002 milliradians (or 3*83 ft) which is very small con- 
sidering also that we started with an initial condition 
value of the X position error of 0.17 milliradians (or 
3256 ft) . 

Figure 69 presents the propagation of Y velocity error. 
It is easily observed that this error drops down to the very 
small value of 0.040 milliradians/hour (or 68 x 10 ft/sec) 
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after one hour and to the negligible error of 0.002 

_c 

milliradians/hour (or 34 x 10 ft/sec) which provides the 
advantage of very accurate evaluation and tracking of the Y 
velocity state variable. 

Similar with the above considerations and even better 
results occur for the case of the X velocity error of the 
integrated I . N . S . /G . P . S . system. We see from Figure 70 that 
this error starting from an initial condition value of 0.3^ 

milliradians/hour (or 2 ft/sec) drops down to 0.1 milli- 

-3 

radians/hour (or 17 x 10 ft/sec) after one hour operation 

and furthermore down to 0.002 milliradians/hour (or 34 x 

-6 

10 ft/sec) after a period of four hours which again 
denotes a very accurate tracking of the X velocity error 
state variable . 

In Figures 71 and 72 the normalized inserted and 
measurement noise of the combined I. N.3./G. P.S. systems are 
presented respectively. 

Thinking of the operation of our system in the long 
term, results of the computer simulation are presented in 
Figures 73 through 79. Using the same input data for our 
I . N . S. /G. P . S. system model and running the program for a 
36 hour process we see that the behavior of the feedback 
configration of the Kalman filter in our system continues to 
be attractive throughout the long term period of interest 
without diverging at any moment. 
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In addition to the above considerations , in order to 
make our system more realistic and compatible to the real 
world’s conditions we put some noise in the two error state 
equations of X and Y position which did not include any 
noise from our theoretical design of the system. So we 
replace the two equations (35-d) and (35-e) in our system of 
equations with the following two equations: 

x(4) = x(6) + AA + K^^CxO) - x(4)] + - x(5)] (35-d') 

and 

x(5) = x(7) + AA + X^^ExCS) - x(4)] + K52Lx(9) - x(5)] (35-e’) 

Assigning to the strength of this intentionally inserted 
noise a value similar to thac of the strength of the gyro 
drift (that is a value of 0.0685 x 10 [rad]“) we ran the 

same program and we achieved results proving that the 
combined I . N . S. /G. P . S. system reacted in a way exactly the 
same as it had reacted without the inserted noise in the X 
and Y position error equations . So we make the conclusion 
that the intentionally inserted noise did not affect the 
operation of our system model neither from the accuracy 
point of view nor from the time point of view. 

The above considerations and results can be seen in 
Figures 80 through 86 for the four (4) hours short term 
process and in Figures 87 through 93 for the 36 hours long 
term operation . 
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In Table VII on the next page we summarized the state 
errors of the combined I. N. S./G. P. S. system after a period 
of two and four hours operation. In the same table we 
included the starting initial condition for each error state 
in order to make our comparisons easier and handy. The 
results included in Table VII are those achieved from the 
computer simulation without any noise corrupting the two 
error states of Y and X position. But since the addition of 
noise with strength similar to that of the gyro drift 
(0.0685 X 10”^ [rad]^) did not affect the system model 
operation as mentioned before, the same Table VII represents 
also the summary of state errors for the real world's system 
model of the I. N . S. /G . P. S. system. 

Up to now we considered our system to be corrupted by 
white Gaussian input noise . Since in the real world in many 
cases the presence of colored noise is apparent we must 
consider the operation of our I . N. S. /G . P. S. system under the 
presence of such noise and compare the results with those 
achieved when the system was driven by white noise. 

In the following section a realistic modeling of the 
I.N.S. component errors is discussed and the results of the 
computer simulation are presented together with the compari- 
son conclusions of the system's operation under colored 
noise versus white noise corruption. 
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TABLE VII 



SUMMARY OF STATE ERRORS FOR THE I . N . S . /G . P . S . SYSTEM MODEL 



State Initial 



Error 


Condition 


Error in 2 hours 


Error in 4 hours 




0.14 mrad 


0.025 mrad 


0.005 mrad 


^E 


0.14 mrad 


0.021 mrad 


0.001 mrad 




0.14 mrad 


0.04 mrad 


0.02 mrad 


SL 


0. 17 mrad 
(3256 ft) 


0.028 mrad 
(536 ft) 


0.01 mrad 
(191.5 ft) 


51 


0.17 mrad 
(3256 ft) 


0.020 mrad 

(383 fft) 


0 . 0002 mrad 
(3.3 ft) 


5il 


0 . 34 mrad/ hr 
(2 ft/sec) 


0 . 040 mrad/hr 
(0.0068 ft/sec) 


0.002 mrad/hr 
(0.00034 ft/sec) 


51 


0 . 34 mrad/ hr 
(2 ft/sec) 


0 . 1 mrad/ hr 
(0.017 ft/sec) 


0.002 mrad/hr 
(0.00034 ft/sec) 
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C. I.N.S. COMPONENT ERROR MODELS 

In Chapter III equation (31) indicates that the I.N.S. 
component errors consist of three gyro drift uncertainties, 
two accelerometer measurement uncertainties, two gyro 
torquer scale factor error's and two geodetic uncertainties. 
Realistic modeling of the two major error components, the 
gyro drift and the accelerometer measurement is described 
below . 

1 . Gyro Drift Uncertainties 



(u)u)jj are each modeled as an exponentially-correlated 
(colored) noise plus an additive random (white) noise: 



where the colored noise is determined by: 




The t^ represents the correlation time of the colored noise 
°i 

and the the strength of the driving white noise obtained 



The three gyro drift uncertainties, (u)u)jj, (u)ujg., 



( u ) oj . = 5 . + w 

, . . °i 



i = N, E, D 



(37) 



using the specified variance 
the formula 



2 



of the colored noise and 



Q, 



V 




( 39 ) 
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2. Accelerometer Measurernent Uncertainties 



The accelerometer measurement uncertaint 
and (u)f£ are modeled in the same way as the gyr 
uncertainties, as colored noise plus white noise 

(u)f^ = a^ + w^ ; i = N, E 

noise of specified 
en by : 



time of the colored 
strength of the driv 

(42) 

Results 

Using the same set of equations (35-a) through 
(35-g) but introducing the appropriate state augmentation 
in order to incorporate the exponentially correlated noise 
for the gyro drift and the accelerometer measurement, we 
simulated the operation of the combined I. N. S./G. P. S. system 
and achieved the following result. 



wnere again w is the wnice 
^i 

and the colored noise is giv 



• 1 

= - r— ''a. 

a . 1 

1 



where z is the correlation 
‘ 2 

with variance a , and the 

^i 

noise v_ is given by: 
a • 

1 



Q = 2 a, /t^ 

V a . a . 

a . II 

1 



3. Computer Simulation 



ies (u)fj^ 
o drift 

(40) 

strength 

(41) 

noise 
ing white 
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The operation of the system proved to be excellent 
for all the used correlation times from 60 seconds up to 
3600 seconds (1 hour). The attitude and navigation errors 
were found to behave in the same way being minimized after a 
period of one hour. Furthermore, the variation of the 



attitude 


and navigation 


errors is 


similar 


with the case of 


the wh 


it 


e noise driven I 


. N. S./G. P 


•S. comb 


ined 


system which 


again 


is 


similar , if not 


ex ac tl y 


the same 


, wi 


th the id 


eal 


I. N. S. 


/G 


. P. S. system . 












In 


Figures 9^ through 100 we 


present 


the 


I. N. S./G. 


P.S. 


system 


0 


peration for an 


ex ponen ti 


ally cor 


rela 


ted input 


noise 


with a 




orrelation time 


of 1 hour 


(3600 5 


ec ) . 


We can 


easily 


see in 


these figures that the beh 


avior of 


the 


combined 




I. N. S. 


/G 


. P. S. system is 


th*e same 


with tha 


t wh 


en white 


noise 


drives 


the input except 


for a ver 


y small 


and 


negligibl 


e 


increa 


se 


of the attitude 


and navigation e 


r ror 


s after the 


f ir st 


ho 


ur of operation . 
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V. CONCLUSIONS AND RECOMMENDATIONS 



A. CONCLUSIONS 

From knowledge gained throughout this work and based on 
the material presented in our study, the following 
conclusions are drawn: 

As far as the I.N.S. errors are concerned we saw chat 

1. The effects of constant gyro drift errors for the 
stationary case are related to Foucault modulation which has 
only a second-order effect on the longitude and latitude 
error states and permit us to neglect it in cheap systems 
designed for navigational purposes. 

2 . For the case of easterly flight, latitude errors 

were reduced by a factor of 1.557 which correponds to the 
ratio Foucault modulation period reduced 

analogously from 33*9 hours to 21.8 hours. 

3. For the westerly flight, case longitude and latitude 
errors grow in approximate proportion to the time-drift rate 
prod uc t . 

4. The accelerometer bias errors have the same effects 
for the stationary case as the gyro drift errors. 

5. The Foucault modulation period increased by the same 
as above factor of 1.557 for the easterly flight case, while 
for the westerly flight decreased by a factor of 0.442 
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corresponding again to the ration for the speciiic 

case . 

6. The error sensitivities remain unchanged for the 
easterly and westerly flight and again we may neglect 
Foucault modulation as producing only second-order effects 
on the navigation states. 

For the combined I. M. S./G. P. S. system the results 
achieved by this study proved that the errors of the 
system’s state variables are damped out in less than one 
hour, denoting effective and successful operation of the 
G.P.S. aiding to the I.N.S. Specifically: 

7 . Using suboptimal Kalman filter gains for one hour 
process, the Y-position error reduced from its initial value 
by a factor of 6 in one hour and by a factor of 17 i-^ four 
ho ur s . 

8. With the same suboptimal Kalman filter gains of one 
hour process, the X-position error proved the system more 
attractive since the error reduced by a factor of 8-5 after 
one hour and by a factor of 856 after four hours. 

9 . Both the X and Y-velocity errors damped out very 
quickly so that after one hour the Y-velocity error reduced 
by a factor of 312.5 and the X-velocity error by a factor of 
117*6, while for a four-hour process both velocity errors 
reduced by a factor of 571 from its initial value. 
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10. The consideration of long term filter's operation 
proved no divergence at all for a process of 36 whole hours. 
The errors remained at the same attractive levels as for the 
four-hour process, fact which enables us to conclude that 
the combined I. N. S. /G. P. S. system works with excellent 
results for both short and long term periods. 





1 1 


. Finally 


the 


operation of the 


comb 


ined 


I. 


N. S./G. P 


. s. 


sys 


tem 


under exponentially correlated 


inpu 


c no 


ise 


proved 


to 


be 


exc 


ellent for 


all 


different correl 


ation 


tiin 


es 


from 60 


sec 


up 


to 


3600 sec , 


wi th 


a negligible inc 


rease 


in 


the 


attitud 


e 


and 


navigation e 


rror 


magnitudes after 


the 


firs 


t hour of 





operation . 

B. RECOMMENDATIONS 

Continued study of this work can be based on the 
following recommendations : 

1. A Kalman filter design study where the primary 
emphasis will be placed upon determination of the "best" 
filter state variable vector. A general covariance analysis 
program for the analysis, evaluation and design of Kalman 
filters, which will help this study, has been tape recorded 
from the Wright Patterson Air Force Base, Air Force Avionics 
Laboratory and modified by the author for use in N.P.S. 
campus computer . 

2. Investigation of various measurement rates using the 
external range measurements from a set of satellites in view 
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among the 18 of the G.P.S. and Kalman filter's performance 
for these rates . 

3 . Possible use of a flight profile generator program, 
which will generate simulated flight patterns instead of 
considering specific only cases for stationary, easterly, 
and westerly flights, together with a satellite motion 
generator required to provide necessary information regard- 
ing the satellites' orbital elements. This recommendation 
applies only to U.S. citizens since such programs already 
exist but they are classified. 

4. Investigation and results evaluation for the effects 
upon filter performance when range-rate measurements are 
available. Then a comparison with the usage of only range 
measurements could be extracted. Another aspect for 
investigation could be the satellite bearing measurements to 
declare best observable satellites and to provide better 
accur ac y . 

5 . Finally a comparison of sequential versus simul- 
taneous measurement would be another area of interest. The 
performance of a filter working with sequential measurements 
is of interest primarily, because of the increased cost of 
equipment required to perform simultaneous measurements and 
computations as compared to the sequential ones. 
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X-SCflLE= 1 . OOE + 01 UNITS INCH, [hours] 

T-SCRLE = 5. OOE-05 UNITS INCH, [rad/meru] 

KWSTflS 

RUN 1 E’N’ vs TIME 

Figure 4. Stationary Case. North Level Error [rad/meru] 
for Constant North Gyro Drift [1 meru] . 
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X-SCflLE= 1 . OOE + 01 UNITS INCH, (hours) 

T-SCRLE = 5. OOE-05 UNITS INCH. [Rad/meru] 

KWSTfiS 

RUN 1 E’E’ VS TIME 



Figure 5. Stationary Case. East Level Error [Rad/meru] for 
Constant North Gyro Drift [1 meru] . 
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X-SCRLE=1 . OOE + 01 UNITS 
T-SCflLE=5. OOE-04 UNITS 



KNSTHS 
RUN 1 



INCH. [hours] 

INCH. [Rad/meru] 

E ’ D ’ VS TIME 



Figure 6. Stationary Case. Azimuth Level Error [Rad/meru] 
for Constant North Gyro Drift [1 meru] . 
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X-3CRLE = 1 . OOE + 01 UNITS INCH, [hours] 
T-SCflLE = 5. OOE-04 UNITS INCH. [Rad/meru] 

KWSTflS 

RUNl DLRVSTIME 

Figure 7. Stationary Case. Latitude Error [Rad/meru] for 
Constant North Gyro Drift [1 meru] . 



87 




X-SCflLE= 1 . OOE + 01 UNITS INCH, [hours] 
T-SCHLE = 2. OOE-03 UNITS INCH. [Rad/meru] 

K W S T fi S 

RUN 2 DLO VS TIME 



Figure 8. Stationary Case. Longitude Error [Rad/meru] for 
Constant North Gyro Drift [1 meru] . 



88 















t 




X-SCRLE= 1 . OOE + 0 I UNITS INCH. Ihours] 

T-SCRLE = 2. OOE-04 UNITS INCH. IRad/hour-meru] 

KWSTflS 

RUN 2 DLfiD VS TIME 

Figure 9. Stationary Case. Latitude Rate Error 

[Rad/hours -meru] for Constant North Gyro Drift 
[1 meru] . 
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X-SCfiLE= ! . OOE + 01 UNITS INCH, [hours] 

T-SCRLE = 5 . OOE-OH UNITS INCH. [Rad/hour*meru] 

KWSTflS 

RUN 2 DLOD VS TIME 

Figure 10. Stationary Case. Longitude Rate Error 

° [Rad/hour ‘ineru] for Constant North Gyro Drift 

[1 meru] . 
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X-SCRLE=1 .OOE + 01 UNITS INCH, [hours] 
T-SCRLE=1 . OOE-01 UNITS INCH. [Rad/meru] 

KWSTflS 

RUN 1 E 'N ’ vs TIME 

Figure 11. Easterly Flight at 600 ft/sec. 

North Level Error [Rad/meru] for Constant 
North Gyro Drift [1 meru] . 
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X-SCflLE= 1 . OOE + 01 UNITS INCH, [hours] 
T-SCflLE=l . OOE-01 UNITS INCH. [Rad/meru] 

KkISTflS 

RUN 1 E'E’ VS TIME 

Figure 12. Easterly Flight at 600 ft/sec 

East Level Error [Rad/meru] for Constant North 
Gyro Drift [1 meru] . 
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X-SCRLE = 1 . OOE + 01 UNITS INCH, [hours] 

T-SCfiLE = 2. OOE + 00 UNITS INCH. [Rad/meru] 

Kwsins 

RUN 1 E 'D ’ vs T IME 

Figure 13. Easterly Flight at 600 ft/sec. 

Azimuth Level Error [Rad/meruj for Constant 
North Gyro Drift [1 meru] . 
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X-SCflLE=l .OOE + Oi UNITS INCH, [hours] 
Y-SCRLE= 1 . OOE + 00 UNITS INCH. [Rad/meru] 

KWSTRS 

RUN 1 Din VS TIME 



Figure 14. Easterly Flight at 600 ft/sec. 

Latitude Error [Rad/meru] for Constant North 
Gyro Drift [1 meru] . 
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X-SCflLE= 1 . OOE + 01 UNITS INCH, (hours) 
Y-SCRLE= 1 . OOE + 0 1 UNITS I NCH . [Rad/meru] 

KWSTflS 

RUN 2 DLO VS TIME 

Figure 15. Easterly Flight at 600 £t/sec. 

Longitude Error [Rad/meru] for Constant North 
Gyro Drift [1 meru] . 
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X-SCflLE= 1 . OOE + 0 1 UNITS INCH, [hours) 
r-SCRLE=1.00E + 00 UNITS INCH. [Rad/hour • meru] 

KWSTflS 

RUN 2 DLflD VS TIME 

Figure 16. Easterly Flight at 600 ft/sec. 

Latitude Rate Error [Rad/hour ‘meruj for Constant 
North Gyro Drift [1 meru] . 
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X-SCflLE= 1 . OOE + 0 1 UNITS INCH. [hours] 

Y-SCPLE=1 . OOE + 00 UNITS INCH. [Rad/hour-meru] 

KNSTnS 

RUN 2 DLOD VS TIME 

Figure 17. Easterly Flight at 600 £t/sec. 

Longitude rate error [Rad/hour ‘meruj for Constant 
North Gyro Drift [1 meru] . 
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X-5CRLE=1 . OOE + 01 UNITS INCH, [hours] 
r-SCRLE=1.00E-01 UNITS INCH. [Rad/meru] 

KNSTflS 

RUN 1 E ’ N ’ VS T I ME 

Figure 18. Westerly Flight at 600 ft/sec. 

North Level Error [Rad/meru] for Constant 
North Gyro Drift [1 meru] . 
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X-SCRLE= 1 . OOE + 01 UNITS INCH, [hours] 
T-SCHLE=1 .OOE-01 UNITS INCH. [Rad/meru] 

KWSTnS 

RUN 1 E 'E ’ vs T IME 

Figure 19. Westerly Flight at 600 ft/sec. 

East Level Error [Rad/meru] for Constant North 
Gyro Drift [1 meru] . 
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X-SCRLE= I . OOE + 01 UNITS INCH.Ihours] 

T-SCRLE = 5. OOE + 00 UNITS I NCH . [Rad/meru] 

KWSTnS 

RUN 1 E 'D ' vs TIME 

Figure 20. Westerly Flight at 600 £t/sec. 

Azimuth Level Error [Rad/meru] for Constant 
North Gyro Drift [1 meru] . 
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X-SCflLE = 1 . OOE + 01 UNITS INCH, [hours] 

T-SCflLE = 5.00E’*‘00 UNITS INCH. [Rad/meru] 

K N 5 T fl S 

RUN 1 DLfl VS TIME 

Figure 21. Westerly Fliglit at 600 £t/sec. 

Latitude Error [Rad/meru] for Constant North 
Gyro Drift [1 meru] . 
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X-SCRLE-1 .OOE + 01 UNITS INCH, [hours] 
T-SCflLE-5. OOE + 00 UNITS INCH. [Rad/meru] 

KNSTnS 

RIJN 2 DLO VS TIME 

Figure 22. Westerly Flight at 600 £t/sec. 

Longitude Error [Rad/meru] for Constant North 
Gyro Drift [1 meru] . 
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X-SCRLE=1 .OOE + 01 UNITS INCH, [hours] 
r - 3 C R L E = 5 . 0 0 E - 0 1 UNITS INCH. [Rad/hour*meru] 

KWSTflS 

RUN 2 DLflD VS TIME 

Figure 23. Westerly Flight at 600 ft/sec. 

Latitude Rate Error [Rad/hour ‘me ru] for Constant 
North Gyro Drift [1 meru] . 
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X-SCnLE= 1 . OOE + 01 UNITS INCH, [hours] 
f - S C fl L E = 1 . 0 0 E 00 UNITS INCH. [Rad/hour -meru] 

K W S T fl S 

RUN 2 DLOD VS TIME 

Figure 24. IVesterly Flight at 600 £t/sec. 

Longitude Rate Error [Rad/hour ‘meru] for Constant 
North Gyro Drift [1 meru]. 
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X-SCRLE= 1 . OOE + 0 1 UNITS INCH, [hours] 

T-SCRLE = 5. OOE-07 UNITS I NCH . [Rad/200ug] 

KWSTfiS 

RUN 1 E 'N ’ vs T I ME 

Figure 25. Stationary Case. North. Level Error [Rad/200yg] 
for Constant North Accelerometer Bias [200ug] . 
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X-SCflLE= 1 . OOE+0 1 UNITS 
T-SCfiLE=5. OOE-07 UNITS 



KNSTRS 
RUM 1 



I N L H . [hours ] 

INCH, [Rad/200wg] 

E’E’ VS TIME 



Figure 26. Stationary Case. East Level Error [Rad/200ygJ 
for Constant North Accelerometer Bias [200iig] . 
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X-5CflLE= 1 . OOE + 0 1 UNITS INCH, [hours] 

T-SCRLE = 5. OOE-07 UNITS INCH. [Rad/200yg] 

KNSTflS 

RUN 1 E Ti ’ vs TIME 

Figure 27. Stationary Case. Azimuth Level Error 

[Rad/200vig] for Constant North Accelerometer 
Bias [200vg] . 
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° 000 001 002 003 004 



X-SCflLE=l .OOE + 01 UNITS INCH. Ihours] 
T-SCRLE=5. OOE-07 UNITS INCH. [Rad/200ug] 

KWSTflS 

RUN 1 DLR VS TIME 

Figure 28. Stationary Case. Latitude Error [Rad/200yg] for 
Constant North Accelerometer Bias [200yg] . 
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X-SCRLE= 1 . OOE + 0 1 UNITS INCH, [hours] 
Y-SCRLE= 1 . OOE-OS UNITS INCH. [Rad/200„g] 

KWSTflS 

RUN 2 DLOVSTIME 

Figure 29. Stationary Case. Longitude Error [Rad/200iJg] 
for Constant North Accelerometer Bias [200pgJ. 
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X-SCflLE=1.00E + 01 UNITS INCH, [hours] 

T -SCALE = 2 . OOE -06 UNITS I NC H . [Rad/hour • ZOOyg] 

KNSTRS 

RUN 2 DLflD VS TIME 

Figure 30 . Stationary Case. Latitude Rate Error 

[Rad/hour • ZOOygJ for Constant North Accelerometer 
Bias [ZOOyg] . 
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X-SCflLE= 1 . OOE + 0 1 UNITS 
T-SCRLE=5. OOE-06 UNITS 



Kwsins 

RUN 2 



INCH, [hours] 

INCH. [Rad/hour* 200pg] 



DLOD VS TIME 



Figure 31. Stationary Case. Longitude Rate Error 

[Rad/hour • 200yg] for Constant North Accelerometer 
Bias [200ygJ. 
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X-SCfiLE = l . OOE + 0 1 UNITS INCH, [hours] 

T-SCPLE = 5. OOE-07 UNITS I NCH . [Rad/200ug] 

KNSTnS 

RUN 1 E 'N ’ vs TIME 

Figure 32. Easterly Flight at 600 ft/sec. 

North Level Error [Rad/200ug] for Constant 
North Accelerometer Bias [200ygJ . 
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X-SCRLE=1 . OOE + 01 UNITS INCH, [hours] 

T-SCflLE = 5. OOE-07 UNITS INCH. [Rad/200gg] 

Kwsins 

RUN 1 E 'E ’ vs TIME ' 

Figure 33. Easterly Flight at 600 ft/sec. 

East Level Error [Rad/200yg] for Constant North 
Accelerometer Bias [200yg] . 
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X-SCniE= 1 . OOE+Ol 
T-SCflLE=5. OOE-07 



KWSTnS 
RUN 1 



UNITS 

UNITS 



INCH, [hours] ■ 

INCH. [Rad/200ygJ i 

E 'D ' VS TIME 



Figure 34. Easterly Flight at 600 ft/sec. 

Azimuth Level Error [Rad/200ygJ for Constant North 
Accelerometer Bias [200yg] . 
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X-5CRLE=1 . QOE+0 1 UNITS INCH. [hours] 

T - SC RLE = 5 . OOE - 0 7 UNITS I NCH . [Rad/200pg] 

KWSTnS 

RUN 1 DLfi VS TIME 

Figure 35. Easterly Flight at 600 ft/sec. 

Latitude Error [Rad/200ygJ for Constant 
North Accelerometer Bias [200ygJ . 
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X-SCflLE= 1 . OOE + 01 UNITS 
Y-SCRLE= 1 . OOE-06 UNITS 



KWSTRS 
RUN 2 



INCH, [hours] 

INCH. [Rad/200iig] 

DLO VS TIME 



Figure 36. Easterly Flight at 600 ft/sec. 

Longitude Error [Rad/ 200]ig] for Constant North 
Accelerometer Bias [200yg] . 
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X-SCflLE=1.00E + 01 UNITS INCH, [hours] 

T - S C flL E = 2 . OOE - 0 6 UNITS I NC H . [Rad/hour* ZOOyg] 

Kwsins 

RUN 2 DLRD VS TIME 



Figure 37- Easterly Flight at 600 ft/sec. 

Latitude Rate Error [Rad/hour • ZOOyg] for Constant 
North Accelerometer Bias [ZOOyg]. 
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X-SCfiLE= 1 . OOE + 01 UNITS INCH, [hours] 

Y-SCRLE = 5. OOE-06 UNITS I NCH . [Rad/hour-zooug] 

KWSTflS 

RUN 2 DLOD VS TIME 

Figure 38. Easterly Flight at 600 ft/sec. 

Longitude Rate Error [Rad/hour • 20 OygJ for Constant 
North Accelerometer Bias [200yg] . 
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X-SCRLE- 1 . OOE + 0 1 UNITS INCH, [hours] 

T-SCflLE = 5. OOE-07 UNITS I NCH . [Rad/200ygj 

KNSTnS 

RUN 1 VS TIME 

Figure 39. Westerly Flight at 600 ft/sec. 

North Level Error [Rad/200yg] for Constant 
North Accelerometer Bias [200yg], 
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X-SCflLE= 1 . OOE + 01 UNITS INCH, [hours] 

Y-SCHLE = 5. OOE-07 UNITS INCH. [Rad/200vg] 

KWSTflS 

RUN 1 E 'E ' V5 T IME 



Figure 40. \festerly Flight at 600 ft/sec. 

East Level Error [Rad/200yg] for Constant North 
Accelerometer Bias [200yg], 
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X-SCflLE= 1 . OOE+01 UNITS 
T-SCflLE=5. OOE-07 UNITS 



Kwsins 

RUN 1 



INCH. thouTS] 

INCH. [Rad/200;ig] 

E 'D ’ vs TIME 



Figure 41, Westerly Flight at 600 £t/sec. 

Azimuth Level Error [Rad/200yg] for Constant 
North Accelerometer Bias [ZOOyg], 
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X-SCRLE=1.00E+01 
Y-SCflLE=5. OOE-07 



KNSTflS 
RUN 1 



UNITS 

UNITS 



INCH, [hours] 

INCH. [Rad/200yg] 

DLfl VS TIME 



Figure 42. 



Westerly Flight at 600 ft/sec. 

Latitude Error [Rad/200jig] for Constant North 
Accelerometer Bias [200yg] . 
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X-SCfiLE= 1 . OOE + 01 UNITS INCH, [hours] 
T-SCflLE= 1 . OOE-06 UNITS INCH. [Rad/200vg] 

KWSTfiS 

RUN2 DLOVSTIME 

Figure 43. Westerly Flight at 600 ft/sec. 

Longitude Error [Rad/ZOOyg] for Constant North 
Accelerometer Bias [200yg] , 
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Figure 44. Westerly Flight at 600 ft/sec. 

Latitude Rate Error [Rad/hour ♦ 200yg] for Constant 
North Accelerometer Bias [200yg] . 
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Figure 45. Westerly Flight at 600 ft/sec. 

Longitude Rate Error [Rad/hour ■ 200pgj for Constant 
North Accelerometer Bias [200ygJ. 
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X-SCRLE= 1 . OOE + 01 UNITS INCH. [hours] 

T-SCPLE = 5. OOE-02 UNITS I NCH . [Rad/140urad] 

KNSTfiS 

RUN 1 E’N’ VS TIME 

Figure 46. Stationary Case. North Level Error 

[Rad/140yradj for Initial North Level Error 
[140vrad]. 



126 




X-SCFILE = 1 . OOE + 0 1 UNITS INCH. (hours] 

T-SCflLE = 5. OOE-02 UNITS I NCH . [Rad/l40«rad] 

KWSTRS 

RUN 1 E 'E ’ VS TIME 

Figure 47. Stationary Case. East Level Error [Rad/140yradJ 
for Initial North Level Error [140yradJ . 
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Figure 48. Stationary Case. North Level Error [Rad/140iirad] 
for Initial East Level Error [140yrad] . 
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Figure 49. Stationary Case. East Level Error [Rad/140y rad] 
for Initial East Level Error [140vrad] . 
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X-SCflLE= 1 . OOE + O 1 UNITS INCH, [hours] 

T-SCflLE = 5. OOE-03 UNITS INCH. [Rad/l40urad] 

KWSTflS 

RUN 5 E’N’ VS TIME 

50* Stationary Cas6. North Levsl Error [Rad/140yrad.j 
for Initial Azimuth Level Error [140yradj . 



130 




I 



X-SCflLE=l . OOE+Ol 
T-SCflLE=2. OOE-03 

KWSTflS 
RUN 5 



UNITS INCH, [hours] 

UNITS INCH. [Rad/140iirad] , 

« 

E ’ E ' VS TIME 

i 



Figure 51. Stationary Case. East Level Error [Rad/140yrad] 
for Initial Azimuth Level Error [140yrad] . 

« 



131 



I 




X-SCflLE= 1 . OOE + 01 UNITS INCH, [hours] 

T -SC ALE = 5 . OOE - 02 UNITS INCH. [Rad/(2 ft/sec)] 

KWSTflS 

RUN 1 E’N’ VS TIME 

Figure 52. Stationary Case. North. Level Error 

[Rad/C2 ft/sec)] for Initial Latitude Rate Error 
[0 . 345mrad/hour = 2 ft/sec]. 
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Figure 53, Stationary Case. Azimuth Level Error 

[Rad/ (2 ft/sec)] for Initial Latitude Rate Error 
[0.345 mrad/hour = 2 ft/sec]. 
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X-SCFILE= 1 . OOE + 01 UNITS INCH, [hours] 

T-SCflLE = 5 . OOE-02 UNITS INCH. [Rad/C2 ft/sec)] 

KWSTflS 

RUN 1 DLflVSTIME 

Figure 54. Stationary Case. Latitude Error 

[Rad/(2 ft/sec)] for Initial Latitude Rate Error 
[0. 345mrad/hour = 2 ft/sec]. 
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X-SCFILE= 1 . OOE + 01 UNITS INCH, [hours] 

T - 5 CfiL E = 5 . OOE -02 UNITS INCH. [Rad/(2 ft/sec)] 

KWSTflS 

RUN2 DLOVSTIME 

Figur© 55. Stationary Case. Longitude Error 

[Rad/(2 ft/sec)] for Initial Latitude Rate Error 
[0. 345mrad/hour = 2 ft/sec]. 
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X-SCflLE= 1 . OOE + 0 1 UNITS INCH, (hours] 

T-SCfiLE = 2. OOE-01 UNITS INCH. [(Rad/hour)/ 
KWSTflS 2 ft/sec)] 

RUM 2 DLflD VS TIME 

Figure 56. Stationary Case. Latitude Rate Error 

[CRad/hour)/(2 ft/sec) J for Initial Latitude 
Rate Error [0.345 mrad/hour = 2 ft/sec]. 
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X-SCflLE= 1 . OOE + 0 I UNITS INCH.[hoursJ 
T-SCRLE = 2.00E-01 UNITS INCH. [(Rad/hour}/ 
KWSTflS 2 ft/hour)] 

RUN 2 DLOD VS TIME 

Figure 57. Stationary Case. Longitude Rate Error 

[(Rad/hour}/ (2 ft/sec)] for Initial Latitude 
Rate Error [0 . 34 5nirad/hour = 2 ft/sec]. 
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Figure 58. Stationary Case. North Level Error 

[Rad/ (2 ft/sec)] for Initial Longitude Rate Error 
[0 . 345mrad/hour = 2 ft/sec]. 
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X-SCRLE= 1 . OOE + 0 1 UNITS INCH, [hours] 

T-SCflLE = 2. OOE-OE UNITS I NCH . [Rad/(2 ft/sec) ] 

Kwsins 

RUN 3 E 'E ' vs TIME 

Figure 59. Stationary Case. East Level Error 

[Rad/ (2 ft/sec)] for Initial Longitude Rate 
Error [0 . 345mrad/hour = 2 ft/sec]. 
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X-SCflLE= 1 . OOE+01 UNITS INCH. Ihoursj 
Y-SCflLE = 2. OOE-02 UNITS INCH. IRad/C2 ft/sec)] 

K H S T n s 

RUN 3 DLfl VS TIME 

Figure 60. Stationary Case. Latitude Error [Rad/ (2 £t/sec)] 
for Initial Longitude Rate Error 
[0 . 34 5mrad/hour = 2 £t/sec] . 
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T-SCflLE=5. OOE-02 
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Figure 61. 



Stationary Case. Longitude error [Rad/C2 ft/sec'll 
for Initial Longitude Rate Error v 
[0 . 345mrad/hour = 2 ft/sec] 
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X-SCRLE= 1 . OOE + 0 1 UNITS INCH, [hours] 
Y~SCflLE=^1.00E~01 UNITS INCH. [ (Rad/hour) / 

KNSTflS (2 ft/sec)] 

RUN 4 DLflD VS TIME 

Figure 62. Stationary Case. Latitude Rate Error 

[(Rad/hour) / (2 ft/sec)] for Initial Longitude 
Rate Error [0 . 345mrad/hour = 2 ft/sec]. 
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Figure 63. Stationary Case. Longitude Rate Error 

[(Rad/hour) / (2 £t/sec)] for Initial Longitude 
Rate Error [0 . 345mrad/hour = 2 ft/sec]. 
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Figure 64. Theoretical I.N.S./G.P.S. , 4-Hour Process 
North Attitude Error [Rad] . 
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Figure 65. 



Theoretical I .N.S./G.P.S. , 4-Hour Process. 
East Attitude Error [Rad] , 
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Figure 66. Theoretical I .N. S. /G . P . S. , 4-Hour Process 
Azimuth Attitude Error [Rad] . 
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Figure 67. Theoretical I.N.S./G.P.S. , 4-Hour Process. 
Y Position Error [Rad] . 
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Figure 68. Theoretical I . N. S. /G. P. S. , 4-Hour Process. 
X Position Error [Rad]. 
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Figure 69. Theoretical I . N. S . /G . P . S . , 4-Hour Process. 
Y Velocity Error [Rad/hour] . 
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Figure 70. Theoretical I . N . S . /G. P. S . , 4-Hour Process. 
X Velocity Error [Rad/hour]. 
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Figure 71. Normalized Input Noise Versus Time. 
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Figure 72. Normalized Measurement Noise versus Time. 
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Figure 73. 



Theoretical I . N. S . /G. P . S . , 36-Hour Process. 
North Attitude Error [Rad] . 
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Figure 74. Theoretical I . N. S . /G. P. S . , 36-Hour Process. 

East Attitude Error [Rad] . 
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Figure 75. 



Theoretical I . N. S . /G. P. S . , 36-Hour Process. 
Azimuth Attitude Error [Rad]. 
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Figure 76. 



Theoretical I . N. S . /G. P. S. , 36-Hour Process. 
Y Position Error [Rad]. 
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Figure 77. Theoretical I . N. S . /G. P. S . , 36-Hour Process. 
X Position Error [Rad] . 
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Figure 73. Theoretical I.N.S./G.P.S. , 36-Hour Process. 
Y Velocity Error [Rad/hour] . 
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Figure 79. Theoretical I . N . S . /G. P. S . , 36-Hour Process. 

X Velocity Error [Rad/hour] . 
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Figure 81. 



Realistic I . N. S . /G. P . S . , 4-Hour Process. 
East Attitude Error [Rad] . 




161 




<-SCRLE=l . GQE-CO UNKS INCH, [hours] 
i'-SCfiLE = 5.0CE OS UNITS INCH. [Rad] 

KwSTnS 

RUN 1 E<D> VS TINE 

Figure 82. Realistic I .N. S ./G. P.S. , 4-Hour Process. 

Azimuth Attitude Error [Rad] . 
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Figure 83. Realistic I . N. S. /G. P . S. , 4-Hour Process. 
Y Position Error [Rad] . 
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Figure 84. Realistic I . N. S . /G. P . S . , 4-Hour Process. 
X Position Error [Rad]. 
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Figure 85. Realistic I . N. S . /G. P . S . , 4-Hour Process. 
Y Velocity Error [Rad/hour] . 
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Figure 86. Realistic I . N. S . /G. P . S . , 4-Hour Process. 
X Velocity Error [Rad/hour] . 
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Figure 87. Realistic I . N. S . /G . P. S • , 36-Hour Process. 
North Attitude Error [Rad] . 
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Figure 88. Realistic I . N. S . /G . P , S . , 36-Hour Process. 
East Attitude Error [Rad]. 
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Figure 89. Realistic I . N. S . / G . P . S . , 36-Hour Process. 
Azimuth Attitude Error [Rad] . 
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Figure 90. Realistic I . N. S . /G. P . S . , 36-Hour Process. 

Y Position Error [Rad] . 
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Figure 91. Realistic I . N . S . /G . P. S . , 36-Hour Process. 
X Position Error [Rad] . 
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Figure 92. Realistic I . N. S . /G . P . S . , 36-Hour Process. 

Y Velocity Error [Rad/hour] . 
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Figure 93. Realistic I.N.S./G. P.S. , 36-Hour Process. 

X Velocity Error [Rad/h.our] . 
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Figure 94. Realistic I . N . S . /G . P . S . , Exponentially Correlated 
Input Noise. 4-Hour Process. North Attitude 
Error [Rad] . 
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Figure 95. Realistic I .N. S . /G. P . S . , Exponentially Correlated 
Input Noise. 4-Hour Process. East Attitude 
Error [Rad] . 
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Figure 96. Realistic I .N.S./G. P.S. , Exponentially Correlated 
Input Noise. 4-Hour Process. Azimuth Attitude 
Error [Rad] . 
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Figure 97. Realistic I . N. S. /G, P . S . , Exponentially Correlated 
Input Noise. 4-Hour Process. Y-Position Error 
[Rad] . 
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Figure 98. Realistic I . N. S . /G. P . S. , Exponentially Correlated 
Input Noise. 4-Hour Process. X-Position Error 
iRad] . 



178 




X-SCFILE=1 . OOE + 00 UNITS INCH, [hours] 
r-SCFILE = 2. OOE-04 UNITS INCH. [Rad/hour] 

Kwsrns 

RUN 2 DinOOT vs TI'^E 

Figure 99. Realistic I . N. S. /G. P. S. , Exponentially Correlated 
Input Noise. 4-Hour Process. Y-Velocity Error 
[Rad/hour] . 
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Figure 100. Realistic I . N . S . /G. P . S . , Exponentially Correlated 
Input Noise. 4-Hour Process. X-Velocity Error 
[Rad/hour] . 
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APPENDIX A 



A SIMPLE EXAMPLE: 

KALMAN FILTER APPLICATION TO A RADAR POSITION 
AIDED INERTIAL NAVIGATION SYSTEM 



I. INTRODUCTION 

The application of a Kalman filter to a simplified radar 
position aided inertial navigation system was investigated 
as a first step of our study. Since the case appears to be 
easy to understand difficult concepts and the results prove 
the design expectations we present hereafter this simple 
case formulated according to the concepts and the outline of 
Ref. 2. 

The I.N.S. system was modeled as white noise driving a 

2 

1/s plant. Radar measurements were assumed to be available 
and were similarly corrupted by white noise. 

The differential equations describing the system and the 
Kalman filter were numerically integrated to yield the 
response for a wide range of input conditions and system 
noise statistics. Particular attention was paid to 
conditions in which the noise statistics employed in the 
filter were different from the statistics of the noise 
actually driving the system dynamics and measurements. 

For all cases considered, including those for which 
intentional mismatches in the noise statistics were 
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introduced, the filter was found to perform in an entirely 
satisfactory manner. This is the filter reliably and quite 
accurately tracked the system's dynamics even at the 
presence of at times rather severe levels of noise. 

In the section to follow, the theoretical development 

of the Kalman filter equations will be presented. This 
development is based on that given in Chapiter 6 of Maybeck 
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conditions, is presented. 

II. I.N.S. AIDED BY POSITION DATA 

For this problem, the model of the I.N.S. is simply a 
double integration of noise-corrupted acceleration infor- 
mation, as depicted in Fig. 101. The noise w is a white 
Gaussian noise of zero mean and variance Kernel 

E[ w( t) w( t+T ) ] = Q6 (t) 

entering at the acceleration level, and it is meant to model 
the errors corrupting the I.N.S. accelerometer outputs 
(accelerometer biases and noise, platform misalignment. 
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Radar Aid Model 



1 * 
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I 

» 





etc.). The noise-corrupted acceleration is integrated once 
to yield I . N. S. -ind icated velocity (v^), and a second time 
to obtain inerti ally-ind icated position (Tj^) • 

Similarly a simple model for the radar or radio 
navigation aid is the true position (r^) corrupted by noise 
u, which is again white Gaussian with zero mean. 

The two error state variables for this case are: 

5r(t) = r . (t) - r^(t) 

^ '' (A-1) 

6v( t) = V^( t) - V^. ( t) 

The measurement to be presented to the filter is the 
difference between the inertially indicated position and 
that measured by the radar or radio navigation aid. 

From Figure 96 we have: 

z(t) = r.(t) - r^(t) 

= [r^(t) + 5r(t)] - [r^(t) - u(t)] = (A-2) 

= srCt) + u(t) 

This is a measurement of the error <Sr(t) corrupted by noise 
u( t) . 

To establish the state dynamics model for the error 
states, first let us consider the total states r^ and v^ . 
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(A-3) 
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The true position velocity and acceleration are related 



by : 



r ^ 



Vu 



r 



(A-4) 



Subtracting (A-4) from (A-3) and using the error state 
definitions of (A-1) we find the desired relation as: 



• 




r* 






r 1 






5r 
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1 




6r 
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5 V 
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0 




5v 
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- ^ 
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- - 



w 



(A-5) 



The measurement model z can be expressed in terms of 
errors as : 



z(t) = Cl 0] 



5r 

5v 



+ u(t) =Hx+u(t) 



Since we would like to design a Kalman filter for this 
situation we need to solve the RICATI equation as below: 



P = FP + PF^ + GQG^ - Ph'^R"^HP 



(A-6) 



where 



P = 



1 1 


^12 


21 


^22 



and 



P - P 
^12 21 
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Figure 102. Kalman Filter Block Diagram 
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E[w( t) w( t+T) ] = Q5(x) 



E[u( t)u(t+T)] = R6 (t) 



Solving for RICATI equation for our case v;e get: 



P P 

^11 12 



P P 

^21 22 



2Pi 2 - (P^i)‘"/R 



'■22 - ’’ll'’l2''‘’ 



?22 - P-| 



Q - (?^2^ 

* 

For the steady state case where p = 0 we get the 



elements of covariance matrix in terms of Q, R and the 
1 / 4 

(Q/R) which represents the natural frequency of the 
system : 



1 1 



12 



2 g1/4 ,3/4 



q 1/2 , 1/2 __ p 



21 



22 



2 g3/4 ,1/4 

The Kalman filter equation is 
X = F£ + Gu + K(t) [z - H£] 
which in terms of error quantities becomes 

+ K [z - 6r] 



— - 












6r 
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1 




6r 


6V 

* 
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0 




6 V 



(A-6a) 



ratio 



(A-7) 



(A-8) 



(A-9) 
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where K = PH^R"' = 





Pl 2 


^21 


^^22 




P,,/R 

P,2/R 



(A-10) 



and using the results of the RICATI equation solution we oan 
write 



K 



1 

1 




1 

C 

3 

C\J 

1 


1 

OJ 

4 




(^n)" _ 



where [rad/sec] (A-11) 



From the above information we can draw the block diagram of 
the Kalman filter for this system as shown in Figure 102. 

The initial transient behavior of the filter gains K-j 
and ^2 depends on , but they are within a few percent of 
their steady state values (independent of P^) after 'j^t = 2 , 
so a prediction of time to reach steady state would be 
approximately 2 /<^^ seconds [Ref. 2]. 

The filter can be put into either feedforward configura- 
tion or feedback configuration. Since for our study we use 
the feedback configuration we present here the outline and 
the results for this configuration. A block diagram of the 
system is presented in Figure 103 as depicted in Maybeck's 
work [Ref. 2]. This block diagram allows us to write the 
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system's equations which we will simulate numerically to 
achieve the system's performance. 

We define the outputs of the I.N.S. system corrected by 
feedback from the filter as follows : 



r(t) = r^(t) - 6r(t) 
v( t) = v^( t) - 6v( t) 



(A-12) 



which will be a very helpful mathematical tool since the 
most straightforward means of generating feedback 
implementations is to write the system and filter equations 
in terms of corrected system states. (For our case 
corrected I.N.S. states.) 

Then using the equations (A-3) and (A-9) together wicn 
equation (A-8) we can write the macrix fcrm of the system 
equations as : 



/\ — 
r^( t) -6r ( t) 




0 1 




r ^ ( t) -6r ( t) 




0* 


[a^( t) +w( t| 


(t)' 


v^(t)-6v(t) 




0 0 




V . ( t) -6v( t) 

1 J 




1 




1 

(\J 

1 



and finally we get the simple form: 



C z( t) -6r ( t) ] 
(A-13) 



r(t) 




0 1 




r(t) 
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1^3^ ( t ) +W ( ^ ) J 


K^Ct) 


[r^( t)-r( t) ] 




= 












-f 






v(t) 




0 0 




_v(t)_ 
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K2(t) 


(A-14) 



In the next section the programming and simulation 
results of the system is presented. 
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Figure 103. Feedback Kalman Filter Configuration 



III. COMPUTER SIMULATION AND RESULTS 

We simulated the system for different input signals, 
different noise characteristics (zero mean Gaussian noise 
with different variances) to see the effect of the filter 
for error estimation. 

From the block diagram of the system in Figure 103 we 
can write the following set of equations which we will use 
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numbers to represent the noise into the system. Typical 
simulation runs used an integration step size of 0.01 
seconds and a total run time of 36 seconds by which point 
the filter had easily achieved steady-state operation in 
most cases. 

Shown on the next page is a run summary representing the 
various conditions that were tested. For each of six cases, 
the covariances of the measurement noise (R) and process 
noise (Q) are indicated. Note that in a number of cases, 
the noise statistics used in the filter were chosen to be 
different from those characterizing the input noise entered 
into the system. This intentional ’’mismatch" was done to 
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investigate filter performance under conditions where the 
true ’’real world" noise statistics are inadequately or 
poorly known . 

In particular it was desired to determine whether any 
instances of filter "divergence" could be observed as a 
result of the mismatch in system noise statistics. It is 
noteworthy that for all conditions tested, the filter 
performed in an entirely satisfactory manner with no 
evidence of divergence. 

It should be noted here that in Table XI the noise 
covariances Q and R actually represent the statistics of 
Discrete Noise entering the system at the integration 
interval At = 0.01 seconds. That is: 

Q = £[w^w,^'^] 

where - kAt 
R 

As it is pointed out by Bryson and Ho [Ref. 13] the 
numerical integration is a discrete approximation to a 
continuous system whose noise processes have spectral 
densities given by 

ECw(t)w(-c+t) ] = Q'6 (t) 

E[u( t)u(T+t)l = R' 5(t) 
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The relation between Q and Q' and R and R' is according 



to Br yson : 



Q' = Q*At 



R' = R*At 



Also shown in Table XI for each case, are the filter 
natural frequency and the steady-state values of the Kalman 
gains . 

A brief discussion will now be given of the results for 
each of the six cases. Detailed plots of the variables of 
interest for each case are attached and will be referred to 
in the subsequent discussion. 

1 . Case I 

For this case we used R = 10,000 and Q = 1 for 
the filter. The actual noise however is mismatched with 
R^ = 100 and = 1 and thus the filter assumes the 
measurement noise of the radar position data considerably 
higher than the case is actually. Shown in the attached 
plots on Figures 104 and 105 is the type of noise actually 
entered into the system using a Gaussian random number 
generator. Also shown are the histories of the Kalman 
filter gains K^ and K 2 versus time. The performance of the 
filter for this case is outstanding as evidenced by the two 
plots for Case I in Figures 108 and I 09 . Here the estimated 
position out of the filter coinsides with the true position 
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denoting that the filter tracks the system extremely well. 
Among the other attached plots, Figure 112 presents the 
noise corrupted radar position in a very imposing way. 

2. Case II 

This case represents one in which the filter and 
external noise are "tuned" so that the same noise statistics 
are employed with R = 100 and Q = 25. Again the Kalman 
gains are plotted indicating the time of steady-state 
condition in Figure 114 and 115. As it is depicted from 
Figures 116, 117 and 118 the Kalman filter rapidly 
"locks-on" to the true position and velocity and accurately 
tracks the system thereafter. 

3. Case III 

In this case the filter and external noise are 

"tuned" with R = 100 and Q = 1. The system's initial 

2 

conditions now include a 10 ft/sec constant acceleration 
and it was desired to see how well the Kalman filter kept up 
with the changing input. Once again the performance of the 
filter in accurately tracking the system is excellent. This 
can be verified looking at Figures 124 and 125 and 126 and 
127 respectively where we can see the coinsidence of the 
true and estimated position and velocity. 

4 . Case IV 

For this run the filter and external noise are again 
"tuned" but with increased statistics of R = 400 and Q = 50. 
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The attached set of plots in figures 13 I through 138 present 
the system and the filter operation proving the accurate and 
satisfactory tracking of the system 
5 . Case V 

Now the filter perceives the radar measurement noise 
to be higher than it actually is. The statistics used for 



this ca 


se were 


R = 400 and Q = 


50 for the filter while 


for 


the ext 


ernal no 


ise we used R^ = 


50 


and r 


50 . The 




reliability and 


the 


accuracy of 


the 


system 


is again dep 


ic ted 


in the 


attached 


plots for Case 


V in 


Figures 


139 through 


146. 


5 . 


Case VI 
















In the 
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that 
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f i 1 ter . 


Here we used 


= 50 


and Q = 


1 0. 


The radar 


meas 


ur ement 


noise R^ = R 


= 


400 was 


assumed 


the 


same. The 


set 


of plots 


in Figures 


147 



througn 154 indicate very good performance of the filter 
despite the intentional mismatch introduced for the system 
driving noise . 

IV. CONCLUSIONS 

Following the development of Reference 2, simplified 
equations characterizing the Kalman filter were derived and 
numerically integrated to yield the filter response to a 
wide range of input conditions and system noise statistics. 
Particular attention was paid to conditions in which the 
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noise statistics employed in the filter differed from the 
statistics of the noise records actually driving the system 
dynamics and measurements. 

For all cases considered including those for which 
intentional mismatches in the noise statistics were 
introduced, the filter was found to perform in an entirely 
satisfactory and reliable manner. By that is meant that the 
filter very accurately tracked the system dynamics even in 
the presence of at times rather sever levels of noise. 

Examination of typical time histories for the variables 
of interest, showed that the filter Kalman gains K^ and 
rapidly settled to their theoretical steady state values 
within a time short compared to the average run time. This 
was accompanied by the filter range and velocity estimates 
rapidly locking on to the true system position and velocity 
and accurately tracking it thereafter . 

It is concluded then that the Kalman filter configura- 
tion discussed here above performed extremely well over the 
range of conditions considered. 
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TABLE VIII 



COMPUTER RUNS SUMMARY 



Filter True 



CASE 


R 

[ft] 2 


Q 

sec 


«t 

[ft] 2 


sec 


rad 


«,)3S 


<■< 2)3 


I 


10,000 


1 


100 


1 


0.100 


0.141 


0.010 


II 


100 


25 


100 


25 


0.707 


1 


0.500 


III^ 


100 


1 


100 


1 


0.316 


0 . 447 


0.100 


IV 


400 


50 


400 


50 


0.595 


0.841 


0.354 


V 


400 


50 


50 


50 


0.595 


0.341 


0.354 


VI 


400 


10 


400 


50 


0.398 


0.562 


0. 158 



Initial Conditions: 



Position = 200 ft 

Velocity r 50 ft/sec 

Acceleration = 0 

P. . = P^^ = 10 P.^ = P 



1 1 



4 

P22 = lo » 



12 



21 



0 



R 

Q 

«t 

Qt 

(‘'l^ss 
(K^) ss 



Natural frequency = [Q/R]^^^ 

Measurement noise covariance used in filter 
Process noise covariance used in filter 
True measurement noise covariance 
True process noise covariance 
Steady state gain ^ 

Steady state gain K 2 = 



Vor this 
acceleration 10 



case system 
ft/ sec"^ , 



assumed 



to have constant 
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Figure 104. Normalized Input Noise Versus Time. 
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Figure 105, Normalized Measurement Noise Versus Time. 
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Figure 106. Case I, Kalman Filter Gain to Velocity. 
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Figure 107. Case I. Kalman filter gain to Acceleration. 
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Figure 108. 



Case I. True Position Versus Time. 
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Figure 109. Case I. Predicted Position Versus Time. 
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Figure 110. Case I. Predicted Velocity Versus Time. 
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Figure 111. Case I. I.N.S. Indicated Position Versus 
Time . 
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Figure 112. 



Case I. Radar Indicated Position Versus Time. 
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Figure 113. Case I. I.N.S, Indicated Velocity Versus Time 
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Figure 114. Case II. Kalman Filter Gain to Velocity. 
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Figure 115. Case II. Kalman Filter Gain to Acceleration. 
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Figure 116. Case II. True Position Versus Time. 
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Figure 117. Case II. Predicted Position Versus Time. 
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Figure 118. Case II. Predicted Velocity Versus Tine 
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Figure 119. Case II. I.N.S. Indicated Position Versus Time 
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Figure 120. Case II. Radar Indicated Position Versus Time. 
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X-SCFILE = 1 .OOE + 01 UNITS INCH, [hours] 

T-SCflLE = 2. OOE + 00 UNITS INCH, [ft/sec] 
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Figure 121. Case ir. i.n.s. Indicated Velocity Versus 



217 



: 




000 001 002^ ^ "003 ■ ^lo7 



X-SCRLE=1 . OOE + 01 UNITS INCH, [hours] 

T-SCflLE = 2. OOE + 01 UNITS INCH, [ft/sec] 

KWSTflS 

RUN 1 K1 VS TIME 



122. Cas6 III. Kslms-ii Filtsr Gd.in to Velocity. 
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Figure 123. Case III. Kalman Filter Gain to Acceleration. 
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Figure 124. Case III. True Position Versus Time. 
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Figure 125. Case III. Predicted Position Versus Time. 
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Figure 126. Case III. True Velocity Versus Time. 
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Figure 127. Case III. Predicted Velocity Versus Time. 
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Figure 128. Case III. I.N.S. Indicated Position Versus 
Time . 
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Figure 129. Case III. Radar Indicated Position Versus 
Time . 
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Figur6 130. Case III. I.N.S. Indicated Velocity Versus 
Time . 
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Figure 131. Case IV. Kalman Filter Gain to Velocity. 
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Figure 132. Case IV. Kalman Filter Gain to Acceleration. 
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Figure 133. Case IV. True Position Versus Time. 
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Figure 134. Case IV. Predicted Position Versus Time. 
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Figure 135. Case IV. Predicted Velocity Versus Time. 
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Figure 136. Case IV. I.N.S. Indicated Position Versus 
Time . 



232 




000 001 002 003 OOH 



X-SCnLE=l . OOE+01 
T-SCfiLE=5. OOE+02 



KNSTRS 
RUN 3 



UNITS 

UNITS 



INCH. [hours] 

INCH. [ft] 

R-RflDflR VS TIME 



Figure 137. Case IV. Radar Indicated Position Versus 
Time. 
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Figure 1S8. Case IV. I.N.S. Indicated Velocity Versus Time. 
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Figure 139. Case V. Kalman Filter Gain to Velocity, 
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Figure 140. Case V. Kalman Filter Gain to Acceleration, 
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Figure 141. Case V. True Position Versus Time. 
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Figure 142. Case V. Predicted Position Versus Time. 



238 




X-SCflLE=l . OOE + 01 UNITS INCH, [hours] 

T-SCRLE = 2. OOE + 01 UNITS INCH, [ft/sec] 

KWSTflS 

RUN 2 V-HflT VS TIME 



Figure 143. Case V. Predicted Velocity Versus Time. 
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Figure 144. Case V. I.N.S. Indicated Position Versus Time. 
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Figure 145. Case V. Radar Indicated Position Versus Time. 
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Figure 146. Case V. I.N.S. Indicated Velocity Versus Time. 
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Figure 147. Case VI. Kalman Filter Gain to Velocity. 
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Figure 148. Case VI. Kalman Filter Gain to Acceleration. 
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Figure 149. Case VI. True Position Versus Time. 
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Figure 150. Case VI. Predicted Position Versus Time. 



246 




X-SCFILE = 1 . OOE + 01 UNITS INCH, [hours] 

Y-SCRLE = 2. OOE + 01 UNITS INCH, [ft/sec] 

KWSTRS 

RUN 2 V-HflT VS TIME 



Figure 151. Case VI. Predicted Velocity- Versus Time. 



247 




X-SCflLE=l . OOE + 01 UNITS INCH, [hours] 
Y-SCRLE=2. OOE+01 UNITS INCH. Ift] 

KWSTfiS 

RUN 3 R-INS VS TIME 



Figure 152.. Case VI. I.N.S. Indicated Position Versus Time. 
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Figure 153. Case VT. Radar Indicated Position Versus Time. 
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Figure 154. Case VI. I.N.S. Indicated Velocity Versus Time. 
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APPENDIX 3 



SATELLITE GEOMETRY, VIEW AND RANGE ERRORS 

The range measurement equation will first be developed. 
Next a simple program has been formulated to verify the 
"observability" and "suitability" of at least four 
satellites at any given time. 

A. RANGE MEASUREMENT EQUATION 

The range measuring process is char acter i zed by a set of 
equations, called the range divergence equations, which are 
generated by the user from a combination of I.N.S. and 
satellite information. This range measurement process 
involves the comparison of a measured value of range against 
a predicted value of range. 

The measured range to a satellite is determined by 
measuring the incremental phase shift between the 
satellite's clock and that of the control station wnich 
supports the user . These clocks were synchronized at an 
earlier time. The computed range on the other hand is 
obtained from satellite ephemeris data and user I.N.S. 
supplied position information. 

The fact is that both the measured and the computed 
range values contain in general errors; so by subtracting 
the computed range value from the measured one, the 
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different will contain only tlie associated errors. This 



difference is called "the range divergence." A Kalman 
filter can be constructed to estimate the errors and improve 
the accuracy of the raw range data if these errors can be 
modeled as the outputs of linear systems driven by white 
Gaussian noise [Ref. 8]. 

1 . Range Divergence 



in order to avoid the rotational inconvenience of using 
superscripts or subscripts to keep track of which satellite 
is being referred to. The results are identical for any one 
of the 13 satellites and therefore very easily extended. 



The case of a single satellite will be considered, 



E 




£ = User-Satellite position vector 

r = Earth-User position vector 
—a 

£g = Earth-Satellite position vector 



Figure 155. Range Vector Definition 
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The range vector of interest is the vector £ from the user 
to the satellite. It is explicitly related to the two 
vectors r^ and r^ which are defined and illustrated in the 
above figure. From the geometry it follows that 

L-- u - u 



where 

r = |r| r l^g-r^l = L'L ~ ^Is"-a^ ’ ^Is"-a^ 0-2) 

The measured range to the satellite, r', is composed of two 
parts 



r' = r + 6r' (B-3) 

where, r is the true range and 6r ’ is the error in the 
measured range to the satellite. The computed range to tne 
satellite, r", is in a similar way written: 

r" = r + 5r" (B-4) 

where, r is again the true range and <5r" is the error in the 
value of the computed range. 

The quantity then, that is being observed, is the 
difference of these two range values and it is the called 
"range divergence," Ar . 

Ar = r' - r" = 5r' - 6r" (B-5) 
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2. Errors in Measured Range 



To model the range measurement error , a knowledge of 
the various error sources which are contained in the 
measurement is required and fitting these error sources with 
empirical data. The model used in our study is a simplified 
version of the one found in [Ref. M] with additional 
information of [Ref. 5]. It is a linear combination of 
three components for each satellite measuremenc corrupted by 
white Gaussian noise (w). Each of the separate components 
is a linear system which is driven by white Gaussian noise. 
The range measurement error is modeled by: 

6r' = ob + c( 6T^ - sTg) + w (B-6) 

where 

Sb = range bias 
c = the speed of the light 
6T^ = user clock phase error 
6T r satellite clock phase error 
w = measurement noise 



The error in the range measurement due to ionosphere delay 
is assumed to be included in the satellite clock phase/range 
error. This is a function of the elevation angle and on the 
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order of 15 feet as a good approximation. The bias term, 

5b, in the above equation accounts for the minor effect of 
both speed of light bias and tropospheric delay uncertain- 
ties in each one of the four satellite range measurements. 

3 . Errors in Computed Range 

The computed satellite position includes error which 
depends on the ephemeris data errors, while the I.N.3. 
errors account for the uncertainty in the user's position. 

So far we have 



Is” = Is ^ ^Is" 



r^" = r^ + 5r_" 
—a —a —a 



The error equation is obtained now since. 



(B-7) 



(3-8) 



(r")^ 



= r' 



• r 



II 



and by taking the differential of both sides we get 



2r”5r” = r'’*5r” + 5r”*r” 



or 

5r" = -pi- (r"‘5r") = (fF £") * (3-9) 

We can easily notice that the quantity, pTr £" » of the right 
hand side is a unit vector from the user to the satellite. 



Ir 



If 







(3-10) 
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and also that 



5r " 





(B-1 1) 



Finally the computed error is written as 



5r" = • (srg" - sr^**) 



(B-12) 



Without any significant accuracy loss, it is assumed thar. 
the earth-satellite range error, 5r, is approximately zero 
for the following logic sequence. Satellite orbital 
parameters are updated by the ground tracking network on a 
periodic basis and relayed to the user along with the range 
data. This ephemeris data is quite accurate and any 
uncertainties in computed satellite range can be accounted 
for by increasing the satellite clock phase error [Ref. 4]. 

Thus the computed range error can be written 



The computation of the above equation requires values for 
the unit vector from the user to the satellite and also 
current values for the north, east and azimuth I.N.S. 
position error states 




(B-13) 




[aN, aE, aD]^ 
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Since we are dealing with a stochastic process simulation 
the root-mean-squared (RMS) values of the covariance of the 
three position errors are used. 

The final form of the range divergence equation is 
obtained now by substitution of equations (B-6) and (B-13) 
into the general form equation (B-5) 

Ar = i „ * or '* -r c6T,, - coT + <Sb + w (3-14) 

Since at least four satellites are required as observables 
to correct for the three components of position and the 
clock phase (or time difference), a minimum of four range 
divergence equations need to be solved simultaneously. 

B. SATELLITE OBSERVABILITY 

In the following development of equations we will drop 
out the double prime (") for notational convenience. A 
given satellite must be in-view by the user in order to 
obtain certain measurements. It is then required that the 
satellite must be aoove some specified minimum angle of the 
user's horizon to get a useful signal. This minimum angle 
depends upon the capabilities of the user-equipment and can 
be characterized as arbitrary. In our study the nominal 
value of this angle was selected of ten degrees. This 
observability criterion together with the suitable selection 
of 18 satellites as the total number of satellites for 
global coverage, insures that regardless of the user's 



257 



f 





position, a sufficient and reasonable number of satellites 
will always be in-view, from which a "best set" of the 
required four satellites may be chosen. 

In the following a method for determining whether or not 
a satellite is observable is presented. 



From the above figure and for the following calculations: 

E . = minimum angle of elevation for useful signal 

min 



I satellite 



User 



E= 

Center 




Figure 156. Observability-Criterion Geometry 




r 



User-Satellite position vector 



= Earth-Satellite position vector 
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= Earth-User position vector 

z earth to navigation transformation matrix 
The quantity r^ is most readily expressed in the navigation 

—a 

frame as 





I 




- 




0 




0 


3 


0 


= 


•J 




R+h 




R 



where 

R = radius of Earth 



h r altitude of user 



and the superscript, n, denotes the frame in wnich the 
vector is coordinated (navigation frame). 

Since the earth-satellite position vector in the earth 
frame, r ” , is derived from the ground track latitude and 
longitude of the satellite in the orbit generator and is 
readily available, the vector £ from the user to the satel- 
lite coordinatized in the navigation frame is written as 



lln 




r 

—a 



n 



z C 



0 

0 

R 



(B-16) 
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The unit vector along is given by 



i - 1 r" • 
ir ■ r - ’ 



r = 



(B-17) 



From the geometry we observe that the azimuth component of 
this unit vector is evaluated as the sine of the elevation 
angle, E, or the cosine of its complementary angle A. That 
is 



ir 




0 

1 



i = cos A 



(3-18) 



So far the observability criterion, if the unit vectors from 
the user to the satellite expressed in the navigation frame 
are computed , becomes 



or 



or 



A < A 



max 



cos A > cos 



(B-19) 



if i„ > cos A__„, the satellite is observable 

1 Q lu a A 

if i < cos A^^„, the satellite is not observable 

y* ^ Cj V ^ 



Since we arbitrarily selected for our study the minimum 
elevation angle of ten degrees the above criterion requires 
the azimuth component of this unit vector to be greater than 
the value of cos 80° = 0.17^- 
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(B-20) 



i >0. 174 



A = 80° 



Total deployment consists of three rings or "constella- 
tions'* of six satellites each. The satellite orbits are 
assumed to be planar and circular; in fact, the orbital 
speed and altitude and thus the orbital period of all satel- 
lites is assumed constant. Furthermore, since global cover- 
age is desired, the satellites on any of the three rings are 
equally spaced; thus, the circular arc between any two adja- 
cent satellites on a ring subtends a central angle of sixty 
degrees (6 x 60° = 360°). The satellite identification code 
used is a two-digit code, the first digit denoting the 
particular constellation-ring (1, 2 or 3) and the second 
digit indicating the particular satellite (1 through 6) 
among the six on the denoted ring. As an example, satellite 
32 is the second satellite on the third ring. 

In order to specify the orientation of any one satellite 
with respect to the Earth-fixed frame, three parameters are 
required. The most convenient parameters are the Euler 
angle, from which the direction cosines or unit vectors to 
the satellite may be determined. This process is explicitly 
described in [Ref. 11] and here we will use directly the 
result for the unit vector from Earth to satellite in Earth 
coordinates . 



— s 



= [a 



13’ 23’ 



^ 33 ^ 



(B-21 ) 
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where 



= sin c sin 5 






(B-22) 



r con n cos ? - sin n sin c cos 5 



and the three Euler angles are C, n, ?. 

The initial conditions are the constant orbital para- 
meters of the satellite orbits are given in the following 
two tables. Note that m refers to the mth satellite on the 
designated n ring; for example, the entry of 2m represents 
the remaining five satellites of the second ring. The 
missing entries from the table, ( — ), are dependent upon the 
initial values of the Euler angles 5 , n and ; which are 
specified. Finally the latitude and longitude values refer 
to the ground track of the satellites. 



TABLE IX 



ORBITAL DESIGN CONSTANTS 



Orbital Period 



12 hours 

63^ (all three rings) 
1 1 , 000 n . miles 



Angle of inclination 



Altitude 



Separation arc-angle 
Ring spacing arc-angle 



60° (6 X 60° = 360°) 

120° (3 X 120° = 360°) 
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TABLE X 



SATELLITE INITIAL CONDITIONS 

Satellite # Initial Latitude Initial Longitude Angle 



1 1 


0° 


0° 




0° 


Itn 


— 


-- 




-60(m-1) 


21 


0° 


120° 


(E) 


0° 


2m 


— 






-60(m- 1 ) 


31 


0^ 


-120° 


(W) 


0° 


3m 


— 






-60(m- 1 ) 


Now the 


sequence of equations 


required to 


compute 


the unit 



vectors for each of the 18 satellites is presented. 

C = 53° 

n = 120°(n-l ) - uj^gt (3-23) 

c = -60(m-1) + Ct 

where 

ra = satellite designator (1 through 6) 

n = ring designator (1, 2 or 3) 

= rotational speed of earth = 15°/^ = 1°/240 sec 

C = orbital speed of the satellite 

2r • n 

C = G /r ® ^ — = 2. 1 n. miles/sec 

e s T 
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Using the Euler angles the components of the unit 
vectors of the satellite-earth range, £ 3 ^ » computed 

[Ref. 10, 11]. 



a ^2 sin c sin 5 

a 22 = -(sin n cos t + cos n sin ; cos 5) (B-24) 

a^^ cos n cos ; - sin n sin ; cos 5 

The ground track latitudes and longitudes are given by 

2 2 

Latitude 9 = arc-tan (a^^/ (.^ 22 ^ * ^^33^ ^ (B-25) 

Longitude a = arc-tan (“323^^33^ (B-26) 

The required components of the unit vector of the user- 
satellite range, , in navigation coordinates may now be 
computed . 



1 

0 

1 




'o' 




1 

1 


0 


e — s 


0 


= 




R 




R 




1 

Q 

( 











n T 

i_ = r /r = [i , i„ , i_ 1 

“•^n “ N 

The observability-criterion requires 



(B-27) 



(B-28) 

(B-29) 
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0. 174 



(B-30) 



> cos 80 



A sample output from the computer program for the observa- 
bility-criterion is included in the next pages listing the 
satellites which are in-view at a particular time instant. 
The output table presents in seven columns the most 
important calculated data. Column 1 contains the satellite 
number according to the previously specified code. Columns 
2, 3 and 4 contain the north, east, and azimuth components 
of the unit vector under consideration. Columns 5 and 6 
give the ground track latitude and longitude for each 
satellite. Finally column 7 gives the observability result 
denoting with "O.K." the satellites which fulfill the 
criterion and with those satellites which do not 

fulfill the criterion. 
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SATELLITE POSITION AND OBSERVABI L I TV CRITERION TIME = 360 
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SATELLITE POSITION AND OBSERVABILITY CRITERION TIME = i200 S 
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SATELLITE POSITION AND OBSERVABILITY CRITERION TIME = 3600 
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